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1.0  INTRODUCTION 


This  report  provides  the  theoretical  background,  the  numeri¬ 
cal  procedures,  a  user's  description  and  a  comparison  with  test 
results  for  a  PC  based  computer  program  for  computing  traction 
force  components  and  torque  in  an  elastohydrodynamically ,  lubri¬ 
cated  contact.  The  program  is  dubbed  McFRIC  because  it  was 
designed  primarily  to  compute  PRICtion,  including  the  effect  of 
MicroContacts.  The  program  in  fact  provides  a  comprehensive  tri¬ 
bological  assessment  of  a  general  lubricated,  concentrated  con¬ 
tact  under  combined  rolling,  spinning  and  sliding.  It  reflects 
the  joint  effect  of  34  input  variables  in  computing  the  following 
descriptive  characteristics  of  the  contacts 

1.  The  contact  ellipse  dimensions  and  area. 

2.  The  elastohydrodynamic  (EHD)  film  thickness  both  at  the 
plateau  and  at  the  constriction  that  forms  at  the  rear 
of  a  lubricated,  concentrated  contact  under  fully 
flooded  (unstarved)  and  isothermal  lubricant  inlet  con¬ 
ditions. 

3.  A  film  thickness  correction  factor  accounting  for  a 
viscosity  decrease  of  the  oil  in  the  contact  inlet  due 
to  shear  heating. 

4.  A  film  thickness  correction  factor  accounting  for  lubri¬ 
cant  starvation  in  the  contact  inlet. 
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The  apportionment  of  the  applied  load  between  the 
asperities  and  the  lubricant  film,  using  the 
Greenwood-Williamson  microcontact  model  with  parameters 
estimated  from  the  ordinary  output  of  a  stylus  profile 
device. 

The  estimated  mean  square  curvature  of  the  composite 
surface  roughness  profile,  the  density  of  surface  sum¬ 
mits  and  the  mean  asperity  tip  radius. 

The  mean  number  of  asperity  contacts  and  the  real  con¬ 
tact  area,  i.e.,  the  total  contact  area  of  the  elasti¬ 
cally  deformed  asperities. 

The  average  number  of  microcontacts  which  have  undergone 
plastic  flow  within  the  contact  area  at  the  computed 
surface  separation. 

An  index  of  surface  fatigue  behavior  based  on  the  number 
and  area  of  plastic  microcontacts. 

The  magnitude  of  the  spin  torque  and  the  magnitude  and 
direction  of  the  tractive  force  transmitted  between  the 
contacting  bodies  by  the  combined  effects  of  (i) 
shearing  of  the  fluid  film  and  (ii)  coulomb  friction 
between  contacting  asperities.  This  computation  may  be 
performed  isothermally  or  thermally.  With  the  thermal 
option  the  heat  generated  by  both  types  of  interfacial 


friction  raises  the  temperature  of  the  fluid  film  and 
the  surfaces  and  alters  the  lubrica  :t  properties.  The 
computed  torque  and  force  components  are  printed  and 
optionally  plotted  as  a  function  of  sliding  velocity  in 
the  rolling  direction. 

11.  The  average  and  maximum  value  of  the  temperature  of  the 
lubricant  film  and  the  surfaces,  when  the  thermal  option 
is  elected. 

12.  The  power  transmitted  by  the  contact  and  the  power 
dissipated  in  friction.  This  is  printed  and  optionally 
plotted  as  a  function  of  sliding  speed. 

Section  2.0  is  a  narrative  outline  of  the  scope  of  the  trac¬ 
tion  prediction  problem  as  it  is  addressed  in  this  report. 

Section  3.0  is  a  summary  of  the  computational  procedures  used 
for  computing  the  thickness  of  the  lubricant  film  separating  the 
contacting  bodies,  accounting  for  the  effects  of  lubricant  film 
starvation  and  inlet  heating.  The  relationships  employed  to  com¬ 
pute  the  dimensions  of  the  contact  ellipse,  the  total  elastic 
approach  of  the  contacting  bodies,  and  the  contact  pressure  are 
summarized. 

Section  4.0  is  a  description  of  the  methodology  employed  to 
calculate  the  load  sharing  between  the  pressurized  lubricant  film 
and  the  elastically  deformed  surface  asperities  which  penetrate 
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it.  The  analysis  yields  the  density  of  microcontacts,  the  mean 
real  and  mean  apparent  pressure  due  to  asperity  deformation,  the 
density  of  plastic  microcontacts,  and  the  mean  real  elastically 
deformed  area  of  contact.  These  quantities  are  computed  as  func¬ 
tions  of  the  plateau  film  thickness  and  the  RMS  profile  height 
and  slope  of  the  two  contacting  bodies  using  the  assumption  that 
the  surface  roughness  processes  are  isotropic  and  that  the 
spectrum  of  the  composite  roughness  is  a  power  function  of  spa¬ 
tial  frequency.  The  computation  of  the  asperity  contribution  to 
the  total  traction  torque  and  force  components  is  developed 
postulating  coulomb  friction  at  the  microcontacts. 

Section  5.0  is  a  description  of  a  nonlinear  Maxwell  rheologi¬ 
cal  model  adopted  for  computing  fluid  traction.  The 
Trachman-Cheng  constitutive  law  is  used  as  the  nonlinear  viscous 
component.  It  is  shown  how  the  model  is  capable  of  accounting 
for  viscous,  elastic  and  plastic  fluid  behavior  as  appropriate. 

The  numerical  scheme  used  to  solve  for  the  components  of  the 
fluid  shear  stress  as  a  function  of  the  position  within  the  con¬ 
tact  ellipse  is  outlined  in  Section  6.0.  The  fluid  viscosity  is 
taken  to  vary  with  pressure  in  accordance  with  Barns'  equation. 
The  fluid's  limiting  shear  strength  is  taken  to  be  proportional 
to  pressure.  Both  of  these  fluid  properties  are  therefore  spa¬ 
tially  variable  over  the  contact  ellipse  via  the  Hertzian  distri¬ 
bution  of  pressure.  The  integration  of  the  shear  stress 


components  over  the  contact  ellipse  to  yield  the  traction  force 
components,  torque  and  power  loss  is  also  described. 

The  iterative  analysis  whereby  the  steady  state  solid  and 
film  mean  plane  temperatures  are  computed  is  described  in  Section 
7.0.  Fluid  shearing  and  asperity  friction  combine  to  form  the 
heat  source.  This  heat  is  dissipated  by  the  mechanisms  of  con¬ 
duction  and  convection.  The  analysis  accounts  for  the  mean 
effect  of  temperature  in  the  film  on  the  lubricant  viscosity  and 
limiting  shear  strength. 

The  organization  and  logic  flow  of  McFRIC  is  given  in  Section 
8.0  along  with  a  detailed  description  of  the  iterative  thermal 
solution.  This  section  contains  instructions  for  installing  the 
program  and  for  the  preparation  of  input  data.  A  sample  program 
input  file  and  the  corresponding  program  output  are  given. 

Section  9.0  contains  a  description  of  the  traction  tests  con¬ 
ducted  at  AFML  including  lubricant  selection,  specimen  and  test 
rig  description  the  choice  of  test  variables  and  the  specific 
test  matrices  used  for  each  lubricant. 

Lubricant  specific  rules  for  the  computation  of  six  physical 
properties  as  a  function  of  temperature  are  described  for  each  of 
the  three  test  fluids.  Computation  of  the  limiting  shear 
strength  of  the  fluid  and  of  the  shear  modulus  of  the  composite 
system  comprising  the  lubricant  and  the  near  surface  layers  of 
the  contacting  bodies,  is  discussed. 


Reduction  of  the  traction  test  data  and  the  compilation  of  a 
traction  data  base  is  described  in  Section  10.0.  These  data  are 
used  as  input  to  McFRIC  in  an  effort  to  predict  the  experimental 
traction  curves.  The  data  showed  that  the  thermal  effect  was 
overpredicted  when  the  limiting  shear  strength  was  taken  to  vary 
inversely  with  absolute  temperature.  With  this  dependence 
removed  the  relatively  minor  thermal  effects  exhibited  by  the 
data  were  well  explained  with  a  thermally  dependent  viscosity. 

Fits  were  found  to  be  acceptably  good  for  all  cases  in  which 
the  Deborah  number  computed  at  the  mean  Hertzian  pressure  via 
Barus'  Law  exceeded  unity.  This  included  all  of  the  data  for  two 
of  the  oils  and  roughly  half  of  the  data  for  the  third.  For  the 
cases  in  which  the  Deborah  number  was  less  than  unity  the  pre¬ 
dicted  curves  approached  their  asymptotic  value  with  increasing 
strain  rate  too  slowly.  An  approach  is  adopted  for  these  cases 
of  using  an  empirical  viscosity  value  to  bring  the  prediction 
into  accord  with  the  experimentally  observed  traction  curve  slo¬ 
pes.  This  approach  though  successful  is  shown  to  be  equivalent 
to  simply  increasing  the  Deborah  number.  It  is  suggested  that 
the  Trachman-Cheng  model  may  need  to  be  supplemented  with  a 
further  lubricant  dependent  parameter  that  governs  the  speed  of 
convergence  to  the  limiting  shear  strength  as  shear  rate 
increases. 

In  the  concluding  part  of  Section  10.0  regression  fits  to  the 
limiting  shear  strength  and  elastic  modulus  values  for  the  three 


oils  are  listed.  Using  these  approximate  equations  one  may  use 
McFRIC  for  conditions  of  load,  speed,  contact  ellipse  aspect 
ratio  and  temperature  for  which  traction  tests  are  not  available. 

Numerous  contributions  were  made  in  the  early  stages  of  this 
work  by  Gail  Hadden  and  Lea  Sheynin.  Robert  Aman  compiled  the 
traction  data  base  and  authored  the  program  MATPROP.  Mark  Ragen 
converted  McFRIC  to  run  on  a  PC  and  developed  the  input  data 
structure.  He  also  provided  the  information  for  program  users 
given  in  Section  8.0.  John  Walsh  and  Monica  Friday  conducted 
most  of  the  McFRIC  runs  and  reruns  and  performed  the  data 
plotting.  Advice  and  many  valuable  suggestions  offered  in  tech¬ 
nical  discussions  with  Dr.  Luc  Houpert  of  SKF  during  the  course 
of  this  work,  are  gratefully  acknowledged. 
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2.0  PROBLEM  STATEMENT 


The  subject  being  addressed  herein  is  the  lubricated  contact 
of  two  moving  elastic  bodies,  focusing  on  the  problem  of  pre¬ 
dicting  the  resultant  forces  and  torque  due  to  the  tangential 
stresses  distributed  over  their  interface. 

The  forces  are  known  as  friction  or  traction  forces;  the 
torque  as  spinning  moment  or  spinning  torque. 

The  bodies  are  assumed  to  be  bounded  by  surfaces  of  revolu¬ 
tion  with  their  separation  in  the  vicinity  of  a  point  of  defor¬ 
mationless  contact,  assumed  to  be  adequately  approximable  as  a 
second  degree  polynomial  in  a  system  of  coordinates  having  its 
origin  at  the  contact  point.  With  this  assumption,  and  in  the 
absence  of  a  lubricant,  the  equations  of  Hertz  apply  for  the 
calculation  of  1)  the  mutual  approach  of  the  bodies  under  a  load 
P  2)  the  dimensions  a  and  b  respectively,  of  the  semimajor  and 
semiminor  axes  of  the  elliptical  interfacial  contact  area  and  3) 
the  elliptically  paraboloid  distribution  of  interfacial  normal 
pressure.  Figure  2.1  shows  the  contact  and  relevant  velocities 
and  dimensions.  The  two  surfaces  are  presumed  to  move  with 
respective  velocities  uj  and  U2  in  a  direction  coincident  with 
one  of  the  principal  axes  of  the  contact  ellipse.  The  average 
velocity  in  that  direction  is  denoted  u  =  (uj  +  U2>/2.  u  is 

customarily  termed  the  rolling  or  longitudinal  velocity.  The 

i 

difference  Au  =  u^  -  U2  is  known  as  the  sliding  or  slipping 
velocity. 

t 

f 

\ 

\ 
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Av  denotes  the  velocity  difference  in  the  direction  orthog¬ 
onal  to  u.  It  is  known  as  the  transverse  sliding  velocity  or  side 
slip.  Finally/  a  relative  rotational  velocity  known  as  spin  may 
act  about  an  axis  perpendicular  to  the  contact  ellipse  and  through 
its  center. 

Figure  2.1(b)  shows  the  distribution  of  interfacial  pressure 
and  its  equation  expressed  in  an  orthogonal  coordinate  system 
established  at  the  contact  center  and  with  the  x-axis  coincident 
with  the  semimajor  axis  of  the  contact  ellipse.  This  pressure 
integrates  over  the  contact  area  to  a  force  equal  to  the  applied 
normal  force  P. 

If  a  viscous  lubricant  is  introduced  and  adheres  to  the  sur¬ 
faces  as  they  approach  the  converging  inlet  to  the  contact  region, 
a  fluid  pressure  builds  and,  for  the  class  of  fluids  typically 
used  as  lubricants,  the  viscosity  increases  and  the  surfaces 
separate  and  deform  to  allow  the  fluid  to  pass  through  the  con¬ 
tact  region  or  'nip'  as  it  is  referred  to  in  some  literature. 

Figure  2.2  shows  the  consequences  of  the  introduction  of  the 
lubricant  on  the  dry  contact  pressure  distribution  and  shape  of 
the  gap  between  the  contacting  bodies.  The  pressure  distribution 
is  extended  forward  in  front  of  the  dry  contact  ellipse.  The  gap 
is  relatively  flat  except  for  a  constriction  near  the  exit  end  of 
the  contact.  A  large  'spike*  in  the  pressure  distribution  pre¬ 
cedes  the  exit  constriction  as  shown. 


The  exact  thickness  and  slope  of  the  lubricant  film  is  deter¬ 
mined  by  the  joint  solution  of  the  Reynold's  equation  of  hydrody¬ 
namics  and  the  deformation  equations  of  contact  elasticity.  The 
film  separating  the  surfaces  is  consequently  known  as  an  elasto- 
hydrodynamic  (EHD)  film. 

The  thickness  of  the  lubricant  film  separating  the  surfaces 
is  small  (typically  0.1-1  urn)  and  is  frequently  of  the  same 
order  as  the  microscale  roughness  of  the  contact  bodies.  The 
bodies  cannot,  therefore,  in  general  be  considered  to  be  fully 
separated  by  the  lubricating  film.  Instead  there  will  be  a  ran¬ 
dom  number  of  microcontacts  taking  place  through  the  lubricant 
film  as  shown  in  Figure  2.3.  The  expected  number  of  such  con¬ 
tacts  will  depend  for  a  given  contact  size  upon  the  ratio  of  the 
mean  film  thickness  to  the  root  mean  square  value  of  the  lubri¬ 
cant  gap  separating  the  surfaces.  This  ratio  is  often  termed 
the  film  parameter  and  is  designated  by  the  symbol  A.  For  A  >  3 
the  expected  number  of  such  contacts  is  negligible.  This  state 
of  affairs  is  known  as  the  full  film  regime.  For  A  <  3  one 
refers  to  a  partial  film  or  mixed  lubrication  regime. 

In  this  work  we  address  both  the  full  film  and  the  partial 
film  regimes. 

We  may  now  state  the  problem  being  considered. 

Given  the  following  variables: 
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1. 


The  size  and  shape  of  the  contacting  bodies. 

2.  The  material  of  the  contacting  bodies. 

3.  The  characteristics  of  the  surface  roughness  of  the  con¬ 
tacting  bodies. 

4.  The  ambient  temperature. 

5.  The  lubricant  properties  as  functions  of  pressure  and 
temperature . 

6.  The  imposed  kinematic  conditions  uj_,  U2,  Av  and  w, 

calculate  the  magnitude  and  direction  of  the  shear  force  x  at 
each  point  within  the  contact  region  considering  the  locally 
variable  conditions  of  relative  velocity,  pressure,  temperature  and 
gap  width  at  each  point.  Sum  (integrate)  these  forces  to  give 
the  orthogonal  components  Fx  and  Fy  of  the  total  traction  force 
and  the  total  traction  or  spin  moment  M. 

In  the  above  statement  of  the  problem  the  kinematics  have 
been  considered  known.  In  implementation,  e.g.  in  an  analysis  of 
bearing  dynamics,  the  converse  problem  is  met  i.e.  it  is  necessary 
to  determine  the  velocities  (e.g.  Au)  for  which  the  resultant 
traction  force  provides  equilibrium.  A  series  of  solutions  are 
therefore  usually  wanted  giving,  e.g.,  the  traction  force  Fx  in 
the  rolling  direction  or  perhaps  the  traction  coefficient  Fx/P  as 
a  function  of  Au  or,  more  often,  as  a  function  of  the  slip-to- 
roll  ratio  Au/u.  Thus,  a  full  solution  of  the  traction  predic¬ 
tion  problem  should  include  the  facility  of  varying  the 
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kinematics  systematically  and  calculating  the  associated  traction 
force  components  and  spinning  moment. 

The  solution  of  the  traction  prediction  problem  involves 
making  the  appropriate  choice  of  a  relationship  between  the 
conditions  at  each  point  (pressure,  relative  velocity,  tem¬ 
perature,  gap  thickness,  etc.)  and  the  shear  stress  and  in 
deducing  the  constants  and  properties  which  are  embodied  in  that 
relationship. 


3.0  FILM  SHAPE,  CONTACT  DIMENSIONS  AND  PRESSURE  DISTRIBUTION 


In  this  section  the  method  adopted  for  the  computation  of  key 
variables  which  affect  the  solution  of  the  traction  prediction 
problem  is  outlined.  These  variables  are  1)  the  thickness  and 
shape  of  the  lubricant  film  that  separates  the  surfaces,  2)  the 
shape  and  dimensions  of  the  interfacial  contact  area  between  the 
surfaces  and  3)  the  magnitude,  shape  and  location  of  the  inter¬ 
facial  pressure  distribution. 

3 . 1  Film  Shape 

As  noted,  the  introduction  of  lubrication  has  the  effect  of 
separating  the  contacting  surfaces  by  a  film  of  virtually 
constant  thickness  over  a  central  plateau  region  while  pro¬ 
viding  a  somewhat  smaller  film  separation  over  a  narrow  constric¬ 
tion  that  forms  at  the  rear  of  the  contact.  As  is  customary  in 
traction  calculations,  the  film  shape  will  be  taken  to  be 
constant  with  thickness  equal  to  the  computed  plateau  thickness. 

Highly  accurate  solutions  for  the  plateau  film  thickness  have 
been  developed  by  Hamrock  and  Dowson  [1]  for  the  fully  flooded 
condition  where  a  copious  lubricant  supply  is  available  and  for 
the  'starved'  case  where  the  lubricant  meniscus  is  at  a  finite 
distance  Xjj,  forward  of  the  contact  center  as  shown  in  Figure 
3.1.  The  value  of  X5  responds  to  the  method  of  lubricant  supply. 
By  specification  of  X5  it  will  be  possible  for  the  user  of  the 
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traction  model  to  account  quantitatively  for  the  effect  of  lubri¬ 
cant  supply  rate  on  the  traction  condition.  Heating  in  the  inlet 
will  increase  the  inlet  temperature  above  ambient/  thus  lowering 
the  viscosity  and  hence  the  EHD  film.  A  multiplicative  reduction 
factor  is  used  to  account  for  this  effect  as  discussed  below: 

3.1.1  Fully  Flooded  Central  Film  Thickness 

The  plateau  or  central  film  thickness  that  develops  in  a 
lubricated  contact  under  flooded  conditions/  i.e.  conditions  of 
copious  lubricant  supply,  is  calculated  according  to  the  formula 
developed  by  Hamrock  and  Dowson  [1], 

hc,f  =  2*69  Rx  v0*67gi0. 53w~0.067 ( l-o. 61e-0 • 73b/a )  (3.1) 

where 

V  =  n0  u/(E'Rx) 

W  =  P/(E'R2X) 

G1  =  E'a 

oq  =  absolute  viscosity  at  ambient  pressure  and  temperature 

u  =  entrainment  velocity,  i.e.,  the  average  surface  velo¬ 
city  in  the  rolling  direction 

E'  =  2  [  ( ( 1-v2  L  )/E !  +  (1-v22)/E2)]“1  (3.2) 


e1'e2  =  Young's  moduli  of  two  contacting  bodies 


V1 '  v2 

Rx 

a 

a,  b 


Poisson's  ratio  for  the  two  contacting  bodies 
effective  radius  in  rolling  direction 
pressure  viscosity  coefficient 

contact  ellipse  semi-axes  in  the  direction  of  and 
transverse  to  the  rolling  direction 


This  formula  was  developed  by  curve  fitting  to  the  results 
obtained  in  full  computer  solutions  of  the  equations  of  elastic¬ 
ity  and  hydrodynamics.  The  results  were  obtained  for  cases  with 
b  >  a,  i.e.  with  the  contact  ellipse  major  axis  transverse  to 
the  rolling  direction  but,  as  stated  by  Hamrock  and  Dowson, 
remain  plausible  for  a  >  b. 


3.1.2  Starved  Central  Film  Thickness 

The  starved  central  film  thickness  hCfS  is  calculated  by 
multiplying  the  fully  flooded  central  thickness  hC/f  by  the 
'starvation'  factor  $3,0.  4>SfC  *s  computed  as 

I"  "10.29 

♦s,c  "  (Xb/a-l)  I  .  xb/a  <  m*  (3.3) 

I  (m*  -  1)  J 

-1.0  ;  Xb/a  >  m* 

m*  =  1  +  3.06  [(Rx/a)2  (hc, f/Rx) ] °« 58 


(3.4) 


Formulas  comparable  to  Eqs.  (3.1)  and  (3.3)  have  been  deve¬ 
loped  for  the  constriction  film  thickness  as  well.  These  are 
also  incorporated  in  McFRIC  and  the  constriction  film  thickness 
is  printed  for  reference. 

3.1.3  Film  Reduction  Due  to  Inlet  Heating 

Convergence  of  the  lubricant  film  in  the  EHD  contact  inlet 
results  in  heating  of  the  inlet  oil  with  a  consequent  loss  in 
viscosity  and  thinning  of  the  lubricant  film.  Murch  and  Wilson 
[2]  have  analyzed  this  problem  and  derived  a  multiplicative  fac¬ 
tor  <Pt  which  may  be  applied  to  the  isothermally  calculated  film 
to  correct  for  this  effect.  The  factor  is  given  by 

Vt  =  1/(1  +  0.108A0-62)  (3.5) 

where 

A  =  4n0  tiu2/k  (3.6) 

and 

k  =  heat  conductivity  of  the  oil 

a  =  temperature  viscosity  coefficient  (°R)~1  defined  through 
Reynold's  viscosity-temperature  relation 

n  =  ng  e“ti{T“To) 


T  =  temperature 


ng  =  viscosity  at  temperature  Tq  and  ambient  pressure 

u  =  entrainment  velocity 

3 . 2  Contact  Ellipse  Dimensions  and  Pressure 

The  calculation  of  elastic  contact  ellipse  dimensions  and  the 
contact  pressure  follows  classical  Hertzian  theory  [3]. 

Figure  3.2  shows  the  assumed  undeflected  forms  and  dimensions 
of  the  bodies  that  comprise  the  type  of  contacts  that  are  being 
considered.  Prior  to  deflecting  under  the  load  P,  the  boundaries 
of  the  bodies  in  the  vicinity  of  contact  are  surfaces  of  revolu¬ 
tion.  The  principal  planes,  i.e.,  the  orthogonal  planes  in  which 
the  radii  of  curvature  are  largest  or  smallest  are  assumed  to 
coincide.  The  principal  radius  for  Body  I  in  Plane  1  is  denoted 
rji .  Correspondingly,  the  principal  radius  for  Body  II  in  Plane 
1  is  denoted  The  principal  radii  in  Plane  2  are  denoted 

rj2  and  rii2'  respectively. 

Under  the  action  of  the  load  P,  the  surfaces  will  deflect  and 
the  region  of  contact  will  expand  from  a  point  to  an  elliptical 
area.  The  axes  of  the  contact  ellipse  will  be  parallel  to  the 
principal  planes.  Whether  the  major  axis  of  the  contact  ellipse 
is  parallel  to  Plane  1  or  to  Plane  2  depends  upon  the  relative 
magnitudes  of  the  principal  radii. 

It  is  possible  to  consider  a  somewhat  more  general  contact 
situation  in  which  the  principal  planes  for  Bodies  I  and  II  make 


Fig.  3.2  PRINCIPAL  PLANES  &  RADII  FOR  CONTACTING  BODIES 


an  arbitrary  angle  with  each  other.  However,  inasmuch  as  the 
most  complete  state-of-the-art  film  thickness  and  fluid  traction 
models  assume  that  the  rolling  direction  is  parallel  to  one  of 
the  axes  of  the  contact  ellipse,  we  limit  consideration  to  the 
geometry  of  Figure  3.2. 

We  define  the  principal  curvatures  p^,  P2  of  a  body  as  the 
reciprocals  of  the  principal  radii  r^,  r2-  The  curvature  asso¬ 
ciated  with  a  given  direction  is  taken  to  be  algebraically  posi¬ 
tive  if  the  body  is  convex  in  that  direction  and  negative  if 
concave.  The  curvature  sum  for  a  contact  is  defined  as: 


Zp  =  pn  +  pI2  +  Pm  +  Pn2 

The  auxiliary  quantity  F(p)  is  defined  as: 

F ( p )  £  (Pjj  +  Pul)  ~  (Pl2  ♦  PII2 ^ 

Zp 

If  F ( p )  >  0.0,  the  major  axis  of  the  contact  ellipse  is 
parallel  to  Plane  2  in  Figure  3.2.  If  F(p)  <  0.0,  it  is  parallel 
to  Plane  1. 


(3.7) 


(3.8) 


The  ratio  <  of  the  major  to  minor  axes  of  the  contact  ellipse 


is  found  by  solving  the  transcendental  equation 

ABS ( F ( p ) )  =  — E^IC)~^F(K) 

(<  -1)  E(k) 


(3.9) 


where  ABS(*)  denotes  the  absolute  value  function  and  F(x)  and 
E(<)  denote  the  complete  elliptical  integrals  of  the  first  and 
second  kinds  defined,  respectively,  by: 


(3.10) 


F  ( < )  =  /  I1-(1-k  2)sin  1/2  d^> 

0 


E  (  <)  =  /  [  1  -  (  1-k“2  )  s  i  n2  ^  3  1/ 2 

0 


(3.11) 


The  major  axis  i.  i  may  be  calculated  as 


=  (S*  E(<-I-  P)l/3 


(3.12) 


E  ii  Zp 


where 


=  applied  load 


+  i-"n2  l-i 

EI  EII 


reduced  elastic  modulus 


EI,II  “  elastic  moduli  of  Bodies  I  and  II 
UI,II  =  Poisson's  ratio 


The  minor  axis  is  computed  as 


l2  =  a/< 


(3.13) 


The  contact  area  A  is 


A  =  it  l\l2 


(3.14) 


The  total  distance  6  through  which  remote  points  within  the 


contacting  bodies  move  under  the  action  of  the  load  P,  is 


6  =  f(k)  •  f  9Up)  (P/mcE'  )2  ]  V3 
2E(<) 


(3.15) 


The  maximum  pressure  og  acting  on  the  contact  ellipse  is 


o n  =  1.5  P/A 


(3.16 
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The  average  pressure  p  is  simply, 

p  =  P/A  (3.17) 


3 . 3  Pressure  Shift 

As  noted,  the  presence  of  the  lubricant  redistributes  the 
interfacial  pressure  from  its  dry  contact,  Hertzian,  shape. 
Tevaarwerk  and  Johnson  [4]  suggested  that  to  a  first  approxima¬ 
tion  this  effect  could  be  modelled  as  a  forward  shift  6X  of  the 
dry  contact  pressure  distribution.  The  dry  contact  center  still 
serves  as  the  axis  of  spin  since  that  is  kinematically  deter¬ 
mined. 


The  displacement  6X  of  the  load  center  due  to  the  redistribu¬ 
tion  of  Hertzian  pressure  by  the  lubricant  film  as  shown  in 
Figure  3.3,  is  calculable  by  the  formula  developed  by  Hamrock  and 
cited  in  Tevaarwerk  and  Johnson  14] . 

6X  =  4.25  a  (gx ) 0 • 022  (g2)-0.35  (b/a)0.91  (3.18) 

where  g^  and  q2  are  the  following  dimensionless  variables: 


As  noted,  the  magnitude  of  the  contact  ellipse  dimensions  and 


pressure  distribution  are  calculated  using  the  Hertz  equations 


(b)  CONTACT  GEOMETRY 


.  3.3  IDEALIZED  EHD  CONTACT  GEOMETRY  AND 
KINEMATICS  WITH  SHIFTED  PRESSURE 
DISTRIBUTION 


Figure  3.3  shows  the  shifted  pressure  distribution  and  inter¬ 
facial  area  over  which  the  shear  stresses  are  appreciable. 

3. 4  Deborah  Number 

As  discussed  in  depth  below  in  Section  5.0,  the  Deborah 
number  plays  a  key  role  as  a  determinant  of  whether  a  fluid  beha¬ 
ves  viscously  or  elastically  as  it  deforms  under  sliding  con¬ 
ditions.  For  an  isoviscous  fluid,  the  Deborah  number  is  defined 
as 

T  =  nu/aG  (3.21) 

where  n,  u  and  a  as  defined  previously  are  the  viscosity,  entrain¬ 
ment  velocity  and  contact  ellipse  semi  axis  in  the  rolling  direc¬ 
tion.  G  is  the  shear  modulus  of  the  elastic  composite  formed  by 
the  fluid  and  the  contacting  bodies. 

Inasmuch  as  the  viscosity  varies  with  pressure  through  the 
contact,  the  Deborah  number  is  a  point  function  of  the  spatial 
coordinates  within  the  contact  region.  To  compute  a  represen¬ 
tative  number  one  could  use  a  value  of  the  viscosity  averaged  by 
means  of  an  appropriate  pressure  viscosity  law  over  the  contact 
pressure  distribution.  As  a  simple  approach  in  computing  the 
Deborah  number,  the  average  viscosity  is  approximated  as  the 
viscosity  at  the  average  pressure  computed  using  the  Barus  rela¬ 
tion.  That  is, 


n 


(3.22) 


=  n0  e  “P 


where  p,  the  average  pressure  is 


P  =  2/3  o0 


(3.23) 


and  ii0  is  the  viscosity  at  ambient  pressure. 

If  the  Barus'  equation  applies  over  the  full  contact,  the 
average  viscosity  will  substantially  exceed  the  viscosity  com¬ 
puted  in  this  manner.  Nevertheless  the  value  thus  computed  is 
indicative,  in  that  it  responds  to  changes  in  pressure  and 
ambient  viscosity.  Quantitative  interpretation  of  T,  i.e., 
classifying  fluid  behavior  based  on  the  value  of  T  is  discussed 
in  Section  10.0. 
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4.0  ASPERITY  CONTACT  TRACTION  MODEL  AND  SURFACE  FATIGUE  INDEX 

As  indicated  in  Section  2.0,  considering  the  partial  EHD 
regime  brings  the  added  complexity  of  asperity  contact  into  the 
traction  prediction  problem.  Bair  and  Winer  [5]  show  curves  of 
traction  force  against  the  film  parameter  A,  in  which  a  six  fold 
increase  of  traction  occurs  as  A  decreases  from  1.0  to  0.3. 
Clearly  the  role  of  asperity  contact  cannot  be  ignored. 

The  approach  taken  within  program  McFRIC  is  to  use  the 
Greenwood-Wi 11 iamson  (GW)  microcontact  model  [6]  to  determine  the 
average  force  acting  at  a  contacting  asperity,  the  number  of 
microcontacts  per  unit  area  and  the  real  contact  area  as  a  func¬ 
tion  of  the  apparent  area.  A  summary  of  the  GW  model  is  given  in 
Section  4.1  below. 

The  three  GW  model  parameters  are  deduced  in  terms  of  the 
RMS  height  and  RMS  slope  using  the  spectral  estimation  approach 
proposed  by  McCool  [7],  This  method  is  summarized  in  Section 
4.3. 

Assuming  that  coulomb  friction  with  a  known  friction  coef¬ 
ficient  acts  at  the  deformed  asperities  the  statistical  expec¬ 
tation  of  the  net  traction  force  components  and  torque  is 
computed  as  described  in  Section  4.4. 


4.1  The  Greenwood-Will iamson  Microcontact  Model 


The  Greenwood-Wi 1 1 iamson  microcontact  model  applies  to  the 
contact  of  a  rough  surface  and  a  smooth  plane.  It  is  based  on 
the  assumptions  that: 

1.  Asperities  have  spherically  capped  summits  of  constant 
radius  R  irrespective  of  their  height. 

2.  Asperities  are  mechanically  independent,  i.e.  the  load 
they  support  depends  on  their  height  and  not  upon  the 
load  supported  by  neighboring  asperities. 

3.  Asperities  deform  elastically  in  accordance  with  the 
Hertzian  relations  between  deflection,  load  and  contact 
area . 

4.  Summit  height  expressed  as  a  deviation  from  the  mean 
plane  of  the  summits  is  a  random  variable  and  follows  a 
gaussian  probability  distribution  with  standard 
deviation  ag. 

Comparisons  of  the  Greenwood-Williamson  model  with  more 
comprehensive  models  that  relax  many  of  its  seemingly  restrictive 
assumptions,  show  that  it  is  nonetheless  astonishingly  accurate 
in  its  predictions  [8].  In  view  of  this  and  its  simplicity  to 
use,  it  has  been  adopted  as  the  microcontact  model  for  use  within 
McFRIC . 


The  GW  model  requires  three  input  parameters: 


E>SUM  :  the  surface  density  of  summits 

os  :  the  standard  deviation  of  the  probability  distribu¬ 
tion  of  summit  heights  and 

R  :  the  deterministic  (non-random)  radius  of  the  spheri¬ 
cal  summit  caps. 

Given  the  values  of  these  parameters,  the  contact  density 
ZNCON ,  the  total  asperity  supported  force  TOTF,  the  total  real 
area  of  contact  per  unit  nominal  area,  TOTA,  and  the  average 
asperity  supported  force  AVF,  are  computed  as: 

ZNCON  =  Dsum  •  F0(d/os)  (4.1) 

TOTA  =  n  Dsum  R  os  F^d/Og)  (4.2) 

TOTF  =  1.333  A0  DSUM  E'  Rl/2  Og3/2  f3/2  (d/og)  (4.3) 

and 

AVF  =  1.333  E'rV2  0g3/2  F3/2(d/os)/F0(d /og)  (4.4) 

The  density  Np  of  plastic  contacts,  i.e.,  of  contacts  which 
have  experienced  some  degree  of  subsurface  plastic  flow  is  given 
by 

Np  =  DSUM  •  F o  (d/os  +  wp*)  (4.5) 

The  elastically  computed  microcontact  area,  Ap,  of  the 
plastically  deformed  contacts,  per  unit  nominal  area  Ag,  is  given 


Ap/Ao  =  Dsum  •  Fi  (d/0s  +  Wp*) 


(4.6) 


where 


[  (1-v2  i)/^  +  (1-v22)/e23 


-1 


(4.7) 


e1'e2  =  Young's  modulus  for  the  rough  and  smooth  surfaces 


vl'v2  =  Poisson's  ratios  for  the  rough  and  smooth  surfaces 


wp*  h  6.4R(Y/E" )2/os 


(-1.8) 


Y  =  Yield  strength  of  rough  surface  in  simple  tension 


F0  (•  ) , F i (• ) , F 2/2 (* )  =  Functions  defined  in  terms  of  the  standard 

gaussian  density  function.  Values  are 
tabled  e.g.  in  [8].  Piecewise  approxima¬ 
tions  to  these  functions  are  employed 
within  McFRIC. 


4 . 2  Fatigue  Index 

If  the  microcontact  area  is  denoted  A,  the  number  of  plastic 
contacts  acting  over  this  area,  ZNPLAS,  is  given  by: 

ZNPLAS  =  Np  •  A  (4.9) 

The  mean  area  APLAS  of  a  plasticized  asperities  is 

APLAS  =  Ap/A0  •  A/ZNCON  (4.10) 

An  index  of  fatigue  behavior  proposed  in  [19],  is  the 
product  of  the  number  and  average  area  of  plastic  contacts  i.e. , 


9  =  ZNPLAS  •  APLAS  (4.11) 

It  is  reasoned  that  the  higher  the  value  of  <P ,  the  greater 
the  opportunity  for  surface  initiated  fatigue  failure.  No  stu¬ 
dies  have  yet  been  performed  to  quantify  the  relation  between 
fatigue  life  and  4>.  It  is  nonetheless  a  reasonable  index  to  use 
to  gauge  relative  performance  of  a  number  of  surface/lubricant 
combinations. 


4 . 3  Relating  the  GW  Parameters  to  Spectral  Moments 

In  work  subsequent  to  the  development  of  Greenwood  and 
Williamson's  model,  it  was  found  that  the  GW  parameters  could  be 
expressed  in  terms  of  3  quantities,  mg,  m2  and  11)4,  known  as 
spectral  moments.  These  values  may  be  determined  from  a  profile 
z  ( x )  as 

mg  =  AVG  (z2 (x) )  (4.12) 

m2  =  AVG  [  (dz ( x )/dx )2 ]  (4.13) 

m4  =  AVG  [(d2z(x)/dx2  )2  ]  (4.14) 

where  z(x)  represents  the  profile  height  deviation  from  the  sur¬ 
face  mean  plane  at  some  position  x,  relative  to  an  arbitrary  ori- 
gin.  mg  is,  of  course,  Rq  and  is  frequently  called  0  in  the 
Tribology  literature.  m2  is  the  mean  square  slope  and  m4  is, 
very  nearly,  the  mean  square  curvature. 

Under  the  assumption  that  z(x)  is  a  gaussian  random  variable 
the  summit  density,  Dgg^,  is  given  by  Longuet-Higgins  [9]  as 


^SUM  =  (1^4/ n»2 )/ 6*/3 


(4.15) 


Bush  et  al  [10]  suggest  that  the  radius  R  may  be  approximated  as 
the  reciprocal  of  the  average  summit  curvature  and  give  the 
expression 

R  =  0.375  (ir/m4)1/2  (4.16) 

Bush  et  al  also  express  the  summit  height  standard  deviation 
°s  as 

os  =  (1-0.8968/a)1/2  irq1/2  (4.17) 

where  a,  termed  the  bandwidth  parameter,  is  defined  ass 

2 

a  =  (m0m4)/m2  (4.18) 

It  is  more  usual  in  Tribology  to  use  the  profile  mean  plane 
as  a  reference  than  the  summit  mean  plane  used  by  the  GW  model. 
Figure  4.1  shows  that  a  smooth  surface  whose  height  above  the  sum¬ 
mit  mean  plane  is  d,  is  situated  at  a  distance  h  =  d  +  "z  above 
the  profile  mean  plane.  The  difference  in  mean  planes  is 
expressible  as  (Bush  et  al  [10] ) 

z  =  4  (mo/iia)  1/2  (4.19) 

Using  Eqs.  (4.17)  and  (4.19)  the  film  parameter  A  i.e.,  the  ratio 

h/Rq  =  h/o  is  found  to  be  linearly  related  to  d/os  as: 

i 

d/os  =  [(h/a)  -  4/(*a)!/2]/(i-0. 8968/a)l/2  (4.20) 


Fig.  4.1  RELATION  BETWEEN  SUMMIT  AND  SURFACE  MEAN  PLANES 
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For  a  specified  h/a  value  and  given  values  of  mg*  m2  and  1x14, 

one  may  use  Eq.  (4.20)  to  find  d/os  and  thereafter#  Eqs.  (4.1)  to 

(4.4)  to  compute  the  microcontact  conditions  at  that  h/a  value. 

''til  ‘ 

'  [lw 

When  both  surfaces  are  rough  an  equivalent  smooth  and  rough  sur- 

face  is  considered  in  which  the  values  of  mg,  m2  and  104  are 

summed  for  the  two  rough  surfaces,  i.e.: 

| 

mG  =  (mg)i  +  ( mo ) 2  (4.21) 

m 

m2  =  (1x12)1  +  (n»2)2  (4.22) 

1x14  =  {1114 )  1  +  (1114)2  (4.23) 

where  (mn)i  denotes  the  n-th  spectral  moment  for  profile  i 

9 

(i=l,2).  The  mn  values  computed  from  Eqs.  (4.21)  to  (4.23)  are 

■ 

referred  to  as  composite  values.  When  surfaces  are  anisotropic 

there  exist  two  orthogonal  directions,  called  principal  direc- 

tions,  along  which  the  profile  value  of  m2  is  a  minimum  and  a 

maximum.  According  to  [11],  the  value  of  m2,  designated  ( m2>e / 

3|W 

for  an  equivalent  isotropic  surface  may  be  constructed  as  the 

-  ■ 

harmonic  mean: 

!».■ 

(m2)e  =  (m2max  •  m2min)*^  (4.24) 

The  m4  values  in  the  same  two  principal  directions  are  com- 

bined  in  the  same  way  to  give  (m4)e.  In  principle  mg  is  inde- 

jH 

pendent  of  tracing  direction.  If  mg  is  measured  in  two 

directions  along  with  m2  and  m4,  the  ordinary  arithmetic  average 

£ 

and  not  the  harmonic  mean  should  be  used  to  combine  the  two 

values  of  mg. 

1 

n 
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4.4  Estimating  1114  from  Stylus  Profile  Equipment  Output 


As  shown  above,  the  GW  model  is  completely  specified  if  we 
know  mg,  m2  and  1114  for  the  composite  and  for  the  equivalent 
isotropic  surface.  These  quantities  may  be  computed  directly 
from  a  profile  as  the  finite  difference  approximations  to  Eqs. 
(4.12)  to  (4.14).  Alternatively,  mg,  m2  and  1114  may  be  computed 
in  terms  of  the  spectrum  s(f)  of  a  profile  if  the  spectrum  is 
known  or  estimated.  Roughly  speaking  the  spectrum  is  a  function 
that  expresses  how  the  various  frequencies  present  in  a  random 
profile  contribute  to  the  overall  profile  variability. 

It  has  been  found  empirically  that  the  spectrum  of  a  rough¬ 
ness  profile  most  always  plots  nearly  linearly  on  a  set  of 
logarithmic  scales  with  perhaps  spikes  at  frequencies  which 
correspond  to  the  spacing  of  grinding  furrows,  chatter  or  other 
nearly  periodic  features.  It  is  postulated,  therefore,  that  the 
spectrum  is  of  the  form  s(f)  ~  f”^  where  k,  termed  the  spectral 
exponent,  is  the  magnitude  of  the  slope  of  a  plot  of  log  s 
against  log  f.  Introducing  a  proportionality  constant  'c'  gives 

s ( f )  =  cf-k  (4.25) 

It  is  further  postulated  that  the  spatial  frequencies  present 
in  the  profile  are  confined  to  a  passband  ranging  from  a  lower 
frequency  f^  to  an  upper  frequency  f2*  The  lower  frequency  is 

associated  with  the  long  wave  cutoff  of  a  profile  instrument. 
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The  upper  frequency  £2  is  determined  by  the  electronic  filter  of 
the  stylus  instrument  or  by  the  finite  stylus  radius,  whichever 
results  in  a  lower  frequency  value. 

The  n-th  spectral  moment  mn  is  given  in  terms  of  the  spectrum 

by : 

f  2 

mn  =  (2*)n  /  fn  •  s ( f )df  (4.26) 

fl 


Using  Eqs.  (4.25)  and  (4.26)  gives  the  results  tabled  below 
for  the  values  corresponding  to  n=0,  2  and  4,  needed  for  the  GW 
microcontact  analysis. 


m0  =  ctn(f 2/f 1 ) ;k=l 

=  [c/( 1-k ) 1 (f 21-k-f 1 1-kl Jk*l 


(4.27) 


m2  =  c(2n  )2Jln(f  2/f  1)  ;k=3 

=  [c(2*)2/(3-k)l [f23-k_fl3-k] .**3 

m4  =  c(2ir)4jm(f2/fi);k  =  5 

=  [c(2n)2/(5-k)] (f 25_k-f l5-kJ ;k*5 


(4.28) 


(4.29) 


t 

! 

t 

l 

I 


The  quantities  Rq  and  Aq  provided  by  current  generation  pro¬ 
file  instruments  are  related  to  mg  and  m2: 


Rq  =  m0l/2  (4.30) 

I 

Aq  =  m21/2  (4.31) 


38 


Thus,  using  Eqs.  (4.27)  and  (4.28)  the  ratio  Rq/Aq  may  be 
expressed  as: 


Rq/Aq  =  1/2* 


( 3 — k ) 


(1-k) 


[f21-k 


(f 23“k 


1/2 

;  k  H  or 


3 


(4.32) 


For  specified  and  f2  and  measured  values  of  Rq  and  Aq,  one 
may  solve  Equation  (4.32)  numerically  for  k.  With  k  thus  deter¬ 
mined,  one  may  then  express  the  ratio  m4/n»Q  from  Eqs.  (4.27)  and 
(4.28)  in  terms  of  f]_,  £2  and  k  and  thereby  estimate  1114  as: 


1114 


(2*)4 


(1-k) 


(5-k) 


(k+1  or  3)  (4.33) 


Having  thus  determined  mg#  m2  and  1:14,  the  GW  parameters  are 
computed  from  Eqs.  (4.15)  to  (4.17). 

In  using  Eq.  (4.32)  to  find  k,  the  composite  Rq  and  Aq  values 
are  formed  following  Eqs.  (4.21)  and  (4.22)  as: 

Rq2  =  Rql2  +  Rq22  (4.34) 

Aq2  =  Aqx2  +  Aq22  (4.35) 

where  the  subscripts  1  and  2  refer  to  the  respective  surfaces. 

The  individual  surface  values  are  input  to  McFRIC.  If  a  sur¬ 
face  is  anisotropic  the  equivalent  Aq  value  should  be  input. 

This  is  computed  from  the  maximum  and  minimum  values  using  Eq. 
(4.24)  as: 
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Aqe  =  [  ( Aq  )maX  *  <Aq>MINJ  1/2 


(4.36) 


4.5  Expected  Values  of  Traction  Force  Components 
and  Torque  at  Asperities _ 

The  contact  ellipse  dimensions  are  assumed  to  have  their  dry 
contact  values.  The  fluid  pressure  at  a  general  position  within 
the  contact  ellipse  is  assumed  to  be  reduced  by  asperity  loading 
but  to  remain  proportional  to  the  dry  contact  pressures 

o'  ( x, y )  =  0oo  [1- (x/a )2  -  (y/b)2]^  (4.37) 

where  a  is  a  constant  of  proportionality. 

The  value  of  a  is  determined  by  requiring  that  the  sum  of  the 
load  supported  by  elastic  asperity  deformation  ( TOT F )  and  the 
fluid  pressure  equilibrate  the  applied  load  P,  i.e. 

(2/3  aaQ  +  TOTF)*ab  =  P  (4.38) 

so  that 

a  =  1.5  (P/irab  -  T0TF]/ao  (4.39) 

The  solution  requires,  of  course,  that  P  >  TOTF  •  nab  and  may  be 
invalid  for  thin  films  and  light  loads. 

For  lubricant  films  that  are  thick  relative  to  the  surface 
roughness  TOTF  =  0  and  9  =  1.0. 


We  assume  further  that  a  Coulomb  friction  law  governs  at  the 
contact  of  a  pair  of  asperities  on  the  two  mating  bodies 
undergoing  macro-EHD  contact.  It  is  not  implied  that  a  Coulomb 
model  is  literally  true;  micro-EHD  effects  may  actually  be  the 
basis  for  the  thin  film  increase  of  traction.  Wedeven  [12]  has 
shown  experimentally  that  micro-EHD  effects  can  be  responsible 
for  a  20%  variation  in  traction.  At  the  present  state  of 
understanding  a  more  refined  model  is  unwarranted. 

Consider  an  asperity  contact  located  at  coordinates  (x,y) 
within  the  contact  ellipse.  The  friction  force  at  this  asperity 
will  act  in  the  direction  of  the  resultant  sliding  velocity  vec¬ 
tor  at  x,y.  The  components  of  the  sliding  velocity  vector  are 

vx  =  Au  -  “ jy  (4.40) 

Vy  =  Av  +  w  (x-6x)  (4.41) 

The  magnitude  of  the  force  vector  is  the  product  of  a  Coulomb 
friction  coefficient,  f,  and  the  normal  force  on  the  asperity 
pair.  The  normal  force  consists  of  the  sum  of  the  integrated 
pressure  at  that  point  in  the  contact  ellipse  and  the  force  to 
deform  the  asperity  pair  elastically  by  the  amount  by  which  they 
interfere  at  the  computed  film  thickness. 

The  force  of  elastic  deformation  is  a  random  variable  having 
a  mean  value  AVF  over  all  contacting  asperities.  The  probability 
of  an  asperity  contact  occurring  within  a  differential  area  dxdy 


within  the  macrocontact  ellipse,  is,  to  a  first  approximation, 
independent  of  coordinate  position. 

This  probability  is  the  product  of  the  number  of  contacts  per 
unit  area  ZNCON  and  the  differential  area  dxdy,  i.e. 

PROB  [Contact  in  dxdy]  =  (ZNCON)  dxdy  (4.42) 

The  expected  normal  force  due  to  elastic  deformation  at  (x,y) 
is  thus 

E ( N )  =  ( AVF  •  ZNCON)  dxdy  (4.43) 

The  expected  friction  force  at  position  (x,y)  is  thus 

E ( F ( x, y )  )  =  f  [AVF  •  ZNCON  +  TOTA  •  0o  (x,y)]  dxdy  (4.44) 

where  a(x,y)  is  the  Hertzian  pressure  distribution.  0  is  the 
fluid  pressure  reduction  factor  and  TOTA  is  the  total  asperity 
area  per  unit  nominal  area. 

The  x  component  of  the  total  expected  friction  force  is  found 
by  integrating  the  x  projection  of  E(F(x,y))  over  the  contact 
ellipse,  i.e. 

E(FX)  =  fft  [AVF*  ZNCON  +  TOTA  •  00  ( x,  y )  ]  f  vx/  ( v2  x  +  vy2  )VS}dxdy 

(4.45) 

The  y  component  is  similarly  defined  but  with  vy  replacing  vx 
in  the  numerator  of  the  integrand. 


The  expected  spinning  torque  due  to  asperity  contact  is 


E  ( T )  =  //x  •  I ( x, y )  vy  dxdy  -  i7y  •  I(x,y)vx  dxdy  (4.46) 

where 

I  ( x,  y )  =  f  [AVF  •  ZNCON  +  TOTA  •  6o  (x,y)]  /  (vx2  +  vy2  )  Vl 

(4.47) 

These  expected  asperity  traction  forces  and  moments  are 
added  to  the  fluid  film  values  computed  as  if  the  surfaces  were 
smooth. 
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5.0  RHEOLOGICAL  MODEL  FOR  FLUID  TRACTION 

I  .  T  - - - - - ~ 

I 

I  5.1  One  Dimensional  Maxwell  Model 

I  - 

A  survey  of  the  literature  on  traction  in  concentrated  con¬ 
tacts  was  made  in  conjunction  with  this  project,  and  a  summary 
account  is  given  in  the  Interim  Report  [13].  A  pivotal  paper  in 
the  traction  literature,  noted  therein,  is  that  of  Johnson  and 
Roberts  [14],  in  which  it  was  ingeniously  established  experimen¬ 
tally  that  a  fluid  in  a  concentrated  contact  can  deform  elasti¬ 
cally  at  sufficiently  high  pressures  and  rolling  velocities.  The 
key  to  this  demonstration  is  the  fact  that  a  lateral  force  deve¬ 
lops  in  the  presence  of  spin  with  an  elastic  fluid,  but  not  if 
the  fluid  deforms  viscously.  This  lateral  force  was  shown  by 
Gentle  and  Boness  [15]  to  be  essential  for  correctly  predicting 
the  steady  state  ball  rotational  axis  in  angular  contact  ball 
bearings. 

To  account  for  the  viscoelastic  response,  Johnson  and  Roberts 
proposed  that  the  lubricant  be  regarded  as  a  Maxwell  fluid  in 
which,  on  the  application  of  a  shear  rate  y,  the  fluid  shear 
stress  x  varies  with  time  t  in  accordance  with  the  differential 
equation: 

( 1/G)  dt/dt  +  x/n  =  r  (5.1) 


Essentially,  this  model  regards  the  strain  rate  f  as  con¬ 
sisting  of  two  components:  an  elastic  strain  rate  Ye  given  by 


Ye  =  1/G  dt/dt 


(5.2) 


and  a  viscous  component  Yv  corresponding  to  ordinary  Newtonian 
viscous  behavior 

Y*v  =  T/n  (5.3) 

I  in  which  y  is  linear  with  t. 

Referring  to  Figure  3.3(a),  wherein  the  lubricant  passes 
through  the  contact  with  velocity  u  =  (u^  +  U2)/2,  the  elastic 
component  of  shear  rate  can  be  written 

Ye  =  (1/G) (dT/dx) (dx/dt)  =  (u/G)  dT/dx  (5.4) 

Using  Eq.  (5.4),  and  assuming  constant  lubricant  properties, 
the  solution  of  Eq.  (5.1)  is, 

t  ( x )  *  hy [l-exp(-G(x+a)/nu)J  (5.5) 

wherein  the  initial  condition  t=0  at  x=-a  has  been  used.  The 
shear  stress  at  the  contact  center,  t(x=0)  approximates  the 
average  shear  stress  over  the  contact.  It  has  the  value, 

t(0)  =  iy [1-exp [-aG/nu] ]  (5.6) 

-  HY  [ 1-exp (-1/T)]  (5.7) 

where  T  =  nu/aG  is  known  as  the  Deborah  number.  For  small  T 

t(0)  -  ny  (5.8) 
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that  is,  the  fluid  behaves  as  a  purely  viscous  Newtonian  fluid. 

For  large  T, 

exp { -1/T )  =  1-1/T  (5.9) 

and  Eq.  (5.7)  reduces  to 

x  ( 0 )  =  ny/T  =  a  G  r/u  (5.10) 

i.e.  the  fluid  acts  purely  elastically. 

For  intermediate  T  values  the  fluid  behaves  viscoelast ically, 
but  in  any  case,  according  to  Eq.  (5.6),  shear  stress  and  hence 
traction  force  increases  linearly  with  shear  rate  Y.  The  Deborah 
number  is  usually  interpreted  as  the  ratio  of  two  time  intervals: 

1)  the  relaxation  time  of  the  fluid,  n/G,  and  2)  the  transit  time, 
i.e.,  the  time  u/a  for  the  fluid  to  travel  a  contact  half  width. 
Large  T  values  represent  the  case  where  the  relaxation  time, 
i.e.,  the  time  to  achieve  the  Newtonian  shear  stress  is  large 
compared  with  travel  time  through  the  contact. 

Johnson  and  Roberts  recognized  that  the  Maxwell  model  of  Eq. 
(5.1)  would  only  apply  at  low  strain  rates  since  the  shape  of 
experimental  traction  curves  invariably  shows  a  nonlinearity  with 
shear  rate. 

5 . 2  Nonlinear  One-Dimensional  Maxwell  Model 

Subsequently,  Johnson  and  Tevaarwerk  [16]  employed  the 
Ree-Eyring  nonlinear  viscous  relation  in  lieu  of  the  linear 


Newtonian  form  within  a  Maxwell  model.  They  thus  deduced  a 
constitutive  relationship  that  was  capable  of  explaining  1)  the 
nonlinear  behavior  of  experimental  traction  curves  taken  at 
various  values  of  the  Deborah  number,  and  2)  the  existence  of  a 
traction  force  transverse  to  the  rolling  direction  in  the  pre¬ 
sence  of  superimposed  spin. 

Johnson  and  Tevaarwerk  noted  that  the  diminishing  rate  of 
increase  of  traction  force  with  sliding  speed  at  high  Deborah 
numbers  "strongly  suggests  an  elastic/plastic  solid".  The  Eyring 
model,  however,  does  not  reach  an  asymptote  with  increased 
sliding,  and  in  [4],  Tevaarwerk  and  Johnson  proposed  instead  the 
use  for  high  Deborah  numbers  of  an  elastic/plastic  model  under 
which  shear  stress  increases  with  sliding  velocity  in  a  linear 
elastic  manner  until  reaching  a  limiting  shear  stress  xc. 

At  low  shear  rates  the  Eyring  based  model  exhibits  Newtonian 
viscous  behavior,  while  the  elastic/plastic  model  exhibits  purely 
elastic  behavior. 

At  high  shear  rates  the  elastic/plastic  model  yields  a 
constant  limiting  shear  stress;  while  the  shear  stress  under  the 
Eyring  model  continues  to  increase  with  increasing  shear  rate. 

In  Section  5.3  below,  the  nonlinear  viscous  model  of  Trachman  and 
Cheng  [17]  is  used  in  the  same  manner  as  Johnson  and  Tevaarwerk 
employ  the  Ree-Eyring  model.  The  result  is  a  single  model 
capable  of  representing  elastic,  plastic,  linear  viscous 
(Newtonian)  and  nonlinear  viscous  fluid  response. 


5.3  The  Trachman-Chenq  Nonlinear  Viscous  Model 


Trachman  and  Cheng  f.17,18]  proposed  a  nonlinear  viscous  model 
which  varies  smoothly  between  Newtonian  behavior  in  which  shear 
stress  and  strain  rate  are  proportional,  to  a  limiting  value  of 
shear  stress  xc  achieved  at  high  strain  rates.  The  model  is: 

Y  =  T/n(l-T/xc)  ;  x  <  xc  (5.11) 

where 

Y  =  strain  rate 
n  =  viscosity 

x  =  shear  stress 
xc  =  limiting  shear  stress 

Expressing  x  in  terms  of  r,  the  model  becomes, 

f  =  nY*/(nY/Tc  +  D  (5.12) 

For  a  Newtonian  fluid 


t  =  HY  (5.13) 

The  T-C  model  thus  gives  a  lower  shear  stress  value  than 
Newtonian  for  all  strain  rates  but  approaches  Newtonian  behavior 
as  Y  approaches  zero. 

Dividing  through  by  nY,  Eq.  (5.12)  becomes 


In  this  form  it  is  seen  that  for  y  large,  t  approaches  tc. 


Thus,  of  itself,  the  T-C  model  accommodates  Newtonian  behavior 
at  low  strain  rates  and  plastic  behavior  at  large  strain  rates. 

For  an  elastic  fluid  as  we  have  seen 

Y*  =  (U/G)  dt/dx  (5.15) 


5 . 4  The  Nonlinear  One-Dimensional  Maxwell  Model 
Using  the  Trachman-Cheng  Viscous  Component 

Following  Tevaarwerk  and  Johnson,  we  express  the  total  strain 
rate  y  as  the  sum  of  an  elastic  strain  rate  Yg  and  a  (nonlinear) 
viscous  strain  rate  Y y: 

Y  =  ye  +  YV  (5.16) 

Using  Eq.  (5.15)  for  Yg  and  the  T-C  model,  Eq.  (5.11),  for  Yy 
leads  after  rearrangement  to  the  following  differential  equation 
for  t 

dt/dx  =  Gy/u  -  Gtc/un  •  t/(tc-t)  (5.17) 

5. 5  Application  to  a  Line  Contact 

In  order  to  evaluate  the  model  implications,  a  line  contact 
(or  a  strip  across  an  elliptical  contact)  of  width  2a  and  with  an 
x  coordinate  established  at  the  contact  center  is  now  considered 
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as  shown  in  Figure  3.3a.  As  with  the  linear  Maxwell  model,  the 
lubricant  properties  are  taken  as  constant  for  this  purpose. 

Separating  the  variables  in  Eq.  (5.17),  and  integrating  gives 
the  relation  between  t  =  t(x)  and  the  coordinate  x. 


t(x)  dT  x 

f  -  =  /  dx  =  x  +  a 

0  [Gy/u  -  Gt  tc/un(Tc-i)]  -a 


(5.18) 


Integration  of  the  left  hand  side  gives,  after  rearrangements 


un  r  1 

- r  £n  - 

riY*/Tc  +  1  )  ^(1-t/tc)-(t/t 


C)/(^Y/TC) 


(5.19) 


uh(t/tc) 

gThy/tT+II 


=  x+a 


Introducing  the  dimensionless  distance  X  *  x/a  and  the 
following  dimensionless  variables: 


T  =  un/aG 


(5.20) 


W  =  nr/T 


(5.21) 


Eq.  (5.19)  becomes 


(T/(W+1)  )An  [1/t (1-t/tc)-(t/tc)/w}]  +  T(t/tc)/(w+1) 
=  X  +  1 


(5.22) 


The  dimensionless  variable  T  is,  of  course,  the  Deborah 


.‘‘yv-v'y ywy  ">“• 


number. 


Eq.  (5.22)  has  been  solved  to  yield  the  dimensionless  shear 
stress  t/c c  as  a  function  of  dimensionless  shear  rate  W  for 
various  values  of  T  and  coordinate  position  X.  Integrating 
across  the  contact,  i.e.  with  respect  to  X,  yields  the  average 
value  r/Tc  as  a  function  of  W  with  T  as  a  parameter.  Figures  5.1 
and  5.2  show  t/tc  as  a  function  of  W  and  of  W/T,  respectively.  It 
is  of  interest  to  show  how  certain  limiting  forms,  viz.  the 
Newtonian  viscous,  the  pure  elastic  and  the  pure  plastic  models 
appear  on  these  grids. 

Newtonian  Viscous  Model 

For  a  Newtonian  purely  viscous  model, 

*  =  ny  (5.23) 

so  that 

T/Tc  =  r,Vtc  =  W  (5.24) 

represents  a  Newtonian  viscous  model.  Since  under  this  model 
t/tc  does  not  depend  on  x,  the  average  value,  t/tc  is  also  aqual 
to  W.  The  value  of  f/tc  f°r  a  Newtonian  viscous  model  thus  plots 
against  W  as  a  line  of  unit  slope.  This  Newtonian  limit  is  shown 
in  Figure  5.1.  On  a  plot  of  t/tc  against  W/T,  a  Newtonian  viscous 
model  will  appear  as  a  straight  line  with  a  slope  equal  to  the 
Deborah  number.  The  curve  for  T  =  0.001  in  Figure  5.2  is  suf¬ 
ficiently  straight  to  suggest  nearly  Newtonian  viscous  behavior. 
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TAU/TAUc 


Fig.  5.2 

t/tc  vs.  dimensionless  strain  rate  w/t 


A  purely  plastic  model  is  represented  by 


or  t/tc  =  1.0  =  t/tc  (5.25) 

This  line  is  shown  on  Figures  5.1  and  5.2. 

Purely  Elastic  Model 

For  a  purely  elastic  model  one  may  integrate  Eq.  (5.14)  to  give 

t  =  GyaX  +  c/u  (5.26) 

where  c  is  a  constant  of  integration. 

Using  the  initial  condition  t  =  0  at  X  =  -1,  gives 

t/tc  =  Gya/uTc  (X+l)  =  (W/THX+l)  (5.27) 

Since  this  expression  is  linear  in  X,  t/tc  is  equal  to  t/tc  at 
the  average  value  of  X,  i.e.  at  X  =  0. 

Thus, 

7/tc  =  W/T  (5.28) 

In  a  plot  of  VTc  against  W/T,  the  pure  elastic  model  will 
thus  have  unit,  slope.  This  pure  elastic  limit  line  is  shown 
drawn  on  Figure  5.2. 
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In  a  plot  against  W,  the  pure  elastic  model  is  a  straight 
line  of  slope  1/T.  The  curve  shown  in  Figure  5.1  for  T  =  50 
appears  straight  to  within  graphical  accuracy  and  thus  represents 
primarily  elastic  behavior. 


Bidimens ional  Constitutive  Law 

Following  Johnson  and  Tevaarwerk  [16],  the  lubricant  film  in 
the  contact  in  which  shear  stresses  are  developed,  is  assumed  to 
be  uniform  and  thin. 


The  flow  due  to  the  pressure  gradient  in  the  EHD  contact  is 
ignored  as  negligible  over  the  high  pressure  region,  with  the 
consequence  (from  the  momentum  equation)  that  the  shear  stress  in 
the  x  and  y  directions  does  not  vary  across  the  thickness  of  the 
film.  In  these  circumstances  the  bidirectional  version  of  the 
constitutive  relation  that  was  adopted  for  the  one  dimensional 
case  in  Section  5.0,  becomes  the  coupled  system  of  differential 
equations : 


udTx 

Gdx 


n(l-XeAc> 


=  y  i 


(5.29) 


— +  XJL  •  Te  =  Yy  (5.30) 

Gdx  xe  n( l-Te/tc) 


where  tx  and  ty  are  the  orthogonal  components  of  the  interfacial 
shear  stress.  te,  the  equivalent  shear  stress,  is  defined  as: 
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■e  =■  <‘x  *  ‘y  >1/2 


(5.31) 


The  orthogonal  components  of  shear  rate,  Yx  and  Yy  are  given 


by 


( Au-wy )/h 


Yy  =  [Av+^)(x-6x)  ]/h 


(5.32) 

(5.33) 


Introducing  the  dimensionless  stresses 

~  ^  /N 

xx  =  Tx/Tc'  Ty  =  Ty/Tc'  Te  =  Xe/Tc 
and  the  dimensionless  strain  rates, 

Wx  =  Yxn/Tc  Wy  =  Yy  n/ T  c 

the  dimensionless  coordinates 
X  =  x/a  Y  =  y/b 

and  the  Deborah  number, 

T  =  un/Ga 

gives  the  transformed  equations, 

A  A  A  A  A 

T(dTx/dX)  +  (TxTe)/Te(l-Te)  =  Wx 

^  A  /s  /\ 

T(dTy/dX)  +  (TyTe)/Te(l-Te)  =  Wy 


(5.34) 


(5.35) 


(5.36) 


(5.37) 

(5.38) 
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6.0  SOLUTION  SCHEME  FOR  CONSTITUTIVE  EQUATIONS 


Figure  6.1  shows  the  contact  ellipse  with  semiaxes  a  and  b 
aligned  with  the  rolling  and  transverse  directions,  respectively. 
The  contact  is  divided  into  strips  parallel  to  the  x  direction 
and  of  thickness  Ay.  Each  strip  in  turn  is  divided  into  rec¬ 
tangular  elements  of  length  Ax.  A  general  strip  with  midline  at 
coordinate  location  y  has  a  half  width  a(y)  given  by: 

a(y)  =  a[l-(y/b)2  ]  V2  (6.1) 

The  pressure  at  coordinate  position  x,y  is 

P  ( x,  y  )  =  oQ  [l-(x/a)2-(y/b)2  ]V2  (6.2) 

where  0O  is  the  maximum  Hertzian  pressure. 

At  a  fixed  slice,  i.e.,  at  a  fixed  y  coordinate,  the  maximum 
pressure  occurs  at  x  =  0,  i.e., 

Pmax(y)  =  °o  U-(y/b2)]VS  (6.3) 

The  distribution  over  that  strip  can  be  written  as 

p(x,y)  =  pmax(y)  •  1 1 —  ( x/ a < y ) ) 2  J  V52  (6.4) 

since  using  (6.1)  and  (6.3)  in  (6.4)  reproduces  (6.2). 

In  solving  the  constitutive  equations  we  focus  on  strips  such 
as  that  shown  in  Figure  6.1  and  use  the  local  dimensionless  x 
coordinate  X  =  x/a(y)  and  the  global  dimensionless  y  coordinate 


Y  =  y/b.  The  equations  are  solved  for  xx  and  ty  along  each 
strip.  These  stresses  are  then  integrated  over  the  strip  to  give 
the  differential  contribution  of  the  strip  to  the  force  component 
and  the  torque  acting  over  the  whole  contact. 

The  solution  scheme  follows  an  ad  hoc  approach  developed  by 
Houpert  [20]  in  conjunction  with  a  different  constitutive  rela- 

A  A 

tionship.  The  derivatives  dxx/dX  and  dTy/dX  are  written  in 
finite  difference  form: 

a  a  a 

dcx/dX  =  (Tx-Tx*)/AX  (6.5) 

dxy/dX  =  (Ty-Ty*)/AX  (6.6) 

where  xx,  iy  are  the  current  values  at  position  X  and  xx*,  Xy* 
are  the  values  at  the  previous  area  element  at  coordinate  posi¬ 
tion  X-AX.  Substituting  into  Eqs.  (5.37)  and  (5.38)  gives,  after 
rearrangement : 


wx  +  T  *x*/AX  =  fx[T/AX  +  1/(1-Te)]  (6.7) 

Wy  +  T  Ty*/AX  =  iy[T/AX  +  1/ ( 1-Te ) ]  (6.8) 

Dividing  (6.7)  by  (6.8)  gives, 


tx/t 

Eq. 


=  a  =  wy  +  T 

r  «  m  A  /i.i 


wx  +  T  tx*/AX 


(6.9) 


6.9)  is  used  to  compute  the  ratio  A  at  each  new  position 


in  terms  of  the  values  at  the  previous  X  position  and  modified  by 
the  values  of  the  strain  rates  Wx  and  Wy  at  the  new  position. 


(6.10) 


Xx  =  Ty  •  A 


Also, 


Te  =  tTx2  +  Ty2  1  ^  =  Tx  1 1+A2  ]  V^2  =  xx.B 


(6.11) 


Using  Eq.  (6.11)  to  eliminate  xx  in  the  original  differential 
equation,  Eq.  (5.37)  gives, 


T  dxe/dX  +  Te/(l-te)  =  BWX 


(6. 12) 


Separating  variables  gives 


Tdxe  =  [BWx-Te/(l-Te)]dX 


(6.13) 


Tdxe/BWx-xe/(l-Te)  =  dX 


(6.14) 


Integrating  over  the  interval  X*  =  X  -  AX  to  X  gives 


J  6  Td"tg/^Wx-xe/ ( 1-Te )]  *  /  dX  =  AX 
■■p*  X-AX 


(6.15) 


The  left  hand  side  integrates  to 


F(te) 


F(Te)-F(xe*) 


(6.16) 


whe  re 


F(te)  =  T/(BWX+1)  fcn  [l-(l+l/BWx)ie]  +  TeT/(BWx+l)  (6.17) 


Given  t  *  at  the  abscissa  value  X*,  the  next  successive  value 


is  found  by  numerically  solving: 


F  ( T_ )  -  F(T  *)  =  AX 


(6.18) 


A  golden  section  search  technique  is  used  to  implement  the  solu¬ 


tion.  Having  thus  found  te,  xe  is  found  as 


Te  =  Tc*Te 


(6.19) 


tx  and  x  are  then  computed  as 


Tx  =  Te/® 


(6.20) 


Tv  =  Tx/A 


(6.21) 


The  solution  scheme  starts  with  xe=0  at  the  inlet  end  of  the 


strip  and  bootstraps  across  the  strip  length  solving  for  the 


stress  on  the  next  successive  element  located  at  the  distance  ax 


further  down  the  strip.  Inasmuch  as  the  lubricant  properties 


depend  on  the  local  pressure  and  temperature  at  each  successive 


location,  they  are  adjusted  to  the  local  conditions  of  each  ele¬ 


ment  using  properties  relationships  discussed  further  below. 
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The  actual  stresses  cx  and  *y  are  obtained  by  multiplying  tx 
and  Ty  by  the  locally  variable  value  of  tc.  The  differential 
torque  contribution  at  location  x,y  is 


DT ( x, y )  =  Tx  *  y  ~  Tv  (X_6x) 


(6.22) 


Tx,ty  and  dT(x,y)  are  integrated  over  the  strip  at  position  y  to 
yield  the  differential  force  components  dFx(y)f  dFy(y)  and  the 
differential  torque  contribution  of  the  strip  dT(y). 


a  (y ) 

dFx (y  )  =  J  xxdx 

-a  (y ) 

a  (y ) 

dFy(y)  =  J  tydx 

-a  (y ) 

a  (y ) 

dT (y )  =  J  dT ( x, y )dx 

-a  (y ) 


(6.23) 


(6.24) 


(6.25) 


Once  these  quantities  are  evaluated  for  each  strip  of  the  contact 
ellipse  they  are  integrated  across  the  strips  to  give  the  total 
fluid  force  and  torque  components: 


b 

=  j  dFx (y ) dy 

-b 


=  J  dFy (y )dy 
-b 


=  J  dT ( y ) dy 
-b 


(6.26) 


(6.27) 


(6.28) 


i-  r*~ 


The  combined  fluid  and  asperity  friction  forces  are  obtained 


■.« 
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by  summing  the  fluid  force  components  and  spinning  torque  from 
Eqs.  (6.26)-(6. 28)  and  the  corresponding  expected  coulomb  fric¬ 
tion  values  from  Eqs.  (4.45)  and  (4.46).  Denote  these  sums  Fx, 
Fy  and  T. 

The  combined  frictional  power  loss  (EL)  may  then  be  computed 


EL  =  Fx  •  Au  +  Fy  •  Av  +  T  •  u» 


(6.29) 


For  some  applications,  most  notably  traction  drives,  the 
power  transmitted  (ET)  through  the  contact  is  of  critical  impor¬ 
tance.  This  is 


ET  -  Fx  •  u 


(6.30) 


The  loss  factor  (LF)  is  computed  as  the  ratio  of  the  power 
loss  to  the  power  transmitted,  i.e. 


L.L  =  EL/ET 


(6.31) 


7.0  HEAT  GENERATION  AND  THERMAL  ANALYSIS 

7 . 1  Heat  Generation 

The  heat  generated  per  unit  volume  is  computed  point-by-point 
along  each  strip  in  the  contact  area.  The  total  heat  input  per 
unit  volume,  Q,  is  the  sum  of  the  contribution  due  to  shearing  of 
the  fluid,  Qf,  and  due  to  the  coulomb  friction  at  the  asperities, 
Qa . 


7.1.1  Fluid  Shear  Heat 

The  fluid  contribution  is  computed  as  the  product  of  the 
resultant  shear  stress 

ie  =  [tx2  +  Ty2]l/2  (7.1) 

and  the  viscous  component  of  the  shear  rate  irv.  The  recoverable 
elastic  shear  is  thus  not  included  in  computing  fluid  heat 
generation.  Using  Eq.  (5.11)  gives 

Yv  =  Te/n(l-xe/Tc)  (7.2) 

so  that 

Of  =  *e  *  *v  (7.3) 

Since  tv  must  not  exceed  y ,  a  numerical  check  is  made.  If 


Y V>Y ,  Of  is  computed  as 


As  noted  in  Section  4.0,  the  quantity  TOTF  computed  using  the 
Greenwood-Williamson  microcontact  model  is  the  asperity  load  per 
unit  area  of  apparent  contact.  With  a  coulomb  friction  coef¬ 
ficient  m ,  the  traction  force  per  unit  area  is  wTOTF. 

Multiplying  by  the  resultant  sliding  velocity  and  dividing  by 
film  thickness  h  gives  the  asperity  generated  heat  per  unit 
apparent  volume,  i.e. 

Qa  =  uTOTF[vx2  +  Vy2 ] 1/2/h  =  M«TOTF»Y  (7.5) 

where  vx  and  vy  are  the  sliding  velocity  components: 

vx  =  Av-u>  y  (4.40) 

Vy  =  Av+w(x-6x)  (4.41) 


7 . 2  Thermal  Analysis 

7.2.1  Temperature  Distribution  in  Solid  and  Film 

It  is  assumed  that  heat  generated  within  the  film  along  a 
strip  in  the  x  direction  as  shown  in  Fig.  6.1  is  conducted  to  the 
two  surfaces  and  carried  by  convection  in  the  direction  of 
rolling.  Conduction  in  the  y  direction  is  neglected  as  negli¬ 
gible.  This  permits  the  independent  analysis  of  the  strips  which 
comprise  the  total  contact. 


It  is  also  assumed  that  both  the  solids  bounding  the  lubri¬ 
cant  film  have  the  same  thermal  properties  and  the  same  ambient 
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temperature  Oo.  The  maximum  film  temperature  thus  occurs  on  and 
is  symmetrical  about  the  center  plane.  The  common  temperature  of 
the  bounding  solid  as(x)  varies  with  coordinate  x  and  may  be 
expressed  as  the  sum  of  the  ambient  temperature  and  the  increase 
ao s(x)  caused  by  heating  in  the  contact.  The  maximum  film  tem¬ 
perature  0c,  in  turn  is  expressed  as  the  sum  of  the  solid  tem¬ 
perature  at  the  given  coordinate  position  and  the  increase  Abc(x) 
of  the  center  film  temperature  over  the  solid  temperature,  i.e.: 


oc(x )  =  o  (x)  +  AO  (x) 


(7.6) 


Oc(x)  =  0Q  +  Abg ( x )  +  Attc(x) 


(7.7) 


Fig.  7.1  shows  the  temperature  profile  of  the  solids  as  a 
function  of  coordinate  position  x  and  the  temperature  profile 
across  the  film  at  a  specific  x  position.  Following  Houpert 
[20,21]  we  assume  that  the  temperature  profile  across  the  film 
has  a  nearly  triangular  shape.  With  this  assumption,  Jaeger's 
equation  [22]  for  computing  A6s(x)  becomes 


aos(x)  =  [^PmCm^m^  1  -1/2  J  -2kA0c/h  [ x— 

-a(y) 


(7.8) 


where 


p  =  density 
c  =  specific  heat 


Fig.  7.1  SURFACE  AND  FILM  TEMPERATURE  PROFILES 


k  =  thermal  conductivity 
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h  =  film  thickness 

a(y)  =  half  width  of  contact  zone  at  coordinate  position  y 

The  subscript  'm'  indicates  that  the  property  is  that  of  the 


solid. 


Houpert's  finite  difference  approximation  for  the  mean  plane 


temperature,  accounting  for  conduction  and  convection  is 


A«c(x)  = 


[  {o/k-Pcu/k  ( Abs-A8*  )/Ax{  h2/ 8  +  EA0*] 


(7.9) 


1  +  E 


where 


E  =  pcuh  /12kAx 


(7.10) 


and  p,c  and  k  are  fluid  properties. 


Like  *x*,  Ab*s  and  A9*c  refer  to  the  previous  interval.  0, 
the  power  dissipated  in  the  film  at  position  x  is  computed  as 
described  in  Sections  7.1.1  and  7.1.2. 


Eqs.  (7.8)  and  (7.9)  must  be  solved  iteratively  until 
mutually  consistent  functions  0s(x)  and  0c(x)  are  found.  A 
further  issue  is  that  the  generated  heat  0  depends  on  the  lubri¬ 
cant  properties  which  themselves  depend  on  the  temperature 
distribution  in  the  film  and,  through  the  pressure,  on  the  posi¬ 
tion  x,y. 
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7.2.2  Temperature  and  Pressure  Dependence  of  Lubricant 


Properties 

The  shear  modulus  G  was  assumed  constant  at  its  average 
value,  independent  of  pressure  and  temperature.  The  lubricant 
viscosity  was  initially  assumed  to  vary  with  temperature  0  and 
pressure  p  as 

n(«,p)  =  nQ  exp [-n(0-eo)  +  ap]  (7.11) 

where  nQ  is  the  viscosity  at  ambient  temperature  0o  and  at 
ambient  pressure,  a  is  the  pressure-viscosity  index  and  t*  is  the 
temperature-viscosity  index.  Eq.  (7.11)  may  be  written  as 

n(«,p)  =  hq  exp  [up]  (7.12) 

where 

ne  -  n0  exp[-ti(b-«0)]  (7.13) 

In  this  form  Eq.  (7.12)  is  known  as  Barus'  equation. 

The  viscosity  averaged  across  the  film  at  position  x  is: 
h/2 

n  =  2/h  /  n  [0  ( z  )  ]  dz  (7.14) 

o 

Considering  a  triangular  temperature  distribution  across  the 
film,  this  yields: 

=  nQ  exp[-tiAtts]  x  [ 1-exp  [-0A0C] ]/dA0c  (7.15) 
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The  critical  shear  stress  was  taken  to  be  proportional  to 
pressure  following  Tevaarwerk  and  Johnson  [4]  and,  following 
Lingard  [26],  to  be  inversely  proportional  to  the  absolute  tem¬ 
perature  of  the  lubricant.  This  may  be  expressed  as: 

Tc  =  ^P/T0  (7.16) 

where  4  is  a  constant  of  proportionality,  and  T0  =  6o+460  is  the 
absolute  ambient  temperature.  At  constant  temperature,  the 
average  value  of  tc  is 

7C  =  4p/T0  (7.17) 

Knowing  tc  at  ambient  temperature  one  may  compute  4  as: 

4  =  Tc/p  •  T0  (7.18) 

Letting  Ts  =  t>s+460  be  the  absolute  temperature  of  the  solid  at 

position  x  and  Tc  =  6c+460  the  corresponding  absolute  temperature 
of  the  center  of  the  film,  the  mean  value  of  tc  assuming  a  linear 
temperature  distribution  across  the  film  is 

Tc  “  WTo/(Tc-Ts)]*n[Tc/Ts]  (7.19) 

This  expression  was  initially  used  within  McFRIC  but  it  was 
found  that  the  thermal  effect  on  traction  was  over  predicted  using 
this  expression.  Taking  tc  to  be  a  function  of  pressure  only 
resulted  in  greatly  improved  fits  (cf.  Section  10.0). 

In  the  numerical  solution  for  temperature,  an  initial  film 
temperature  rise  A6c(x)  was  assumed  by  neglecting  convection  in 
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Eq.  (7.9)  and  taking  the  heat  generation  rate  0  to  be  the  product 
of  ic  and  Au/h. 

This  gives 

ASc(x)  =  TcAuh/8  (7.20) 

and  t c  varies  with  x  because  of  its  dependence  on  p  =  p(x,y). 

This  expression  is  used  in  Eq.  (7.8)  to  yield  A0s(x).  Eq.  (7.8) 
contains  a  singularity  at  V  =  x.  In  the  implementation  the  first 
two  values  of  A0s(x)  on  each  strip  are  taken  to  be  zero. 

For  the  third  point,  the  following  approximate  expression 
developed  by  Houpert  (24]  is  used: 

A«s(x)  =  A[23.OA0c(x-Ax)-5.878Aec(x-2Ax)]k(Ax)V2/3h 
where 

A  a  tnpm  cm  km  ul_1/2 
For  the  fourth  point, 

A0S ( x  )  =  A[O.  5773A0c(x-3Ax)-1.17A0c(x-2Ax)  (7.23) 

+  1 1 A0C ( x-Ax ) ]  x  2k(Ax)2/3h 

For  subsequent  points 

A0s(x)  *  A(I+lOA0c(x-Ax)-4A0c(x-2Ax)2k(Ax)1/2/3h  (7.24) 


(7.21) 


(7.22) 
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where 


x-Ax  -a«c(y) 

I  =  j  —  -  d  v  (7.25) 

-a(y)  /  x - f1 

is  evaluated  numerically  using  Simpson's  rule. 

Expressions  (7.21)  and  (7.22)  are  found  by  integrating  by 
parts  over  the  singularity.  Following  the  computation  of  A0s(x) 
the  lubricant  propertis  are  corrected,  and  the  generated  heat 
recomputed.  Eq.  (7.9)  is  then  used  to  recompute  A8c(x),  now 
allowing  for  convection. 


8.0  PROGRAM  DESCRIPTION,  ORGANIZATION  LOGIC  FLOW  AND  USERS 

INFORMATION 

8. 1  Description  of  Computer  Program 

The  McFRIC  program  computes  the  components  of  sliding  trac¬ 
tion  in  the  rolling  direction  and  the  transverse  direction 
(orthogonal  to  the  rolling  direction)  as  well  as  the  spinning 
torque  for  a  general  elastohydrodynamic  (EHD )  contact  that 
experiences  rolling,  sliding  and  spinning.  It  calculates  the 
power  transmitted  and  the  power  lost  for  each  kinematic  con¬ 
dition.  The  lubricant  shear  force  is  modelled  as 
viscoelastic/plastic  with  temperature  and  pressure  dependent 
viscosity  and  limiting  shear  stress.  Asperity  friction  is 
assumed  to  be  Coulombic  and  is  proportional  to  the  asperity  load 
computed  via  the  Greenwood-Williamson  microcontact  model. 

The  program  solution  can  be  isothermal  or  it  can  consider  the 
generation  of  heat  due  to  fluid  shear  and  asperity  friction  and 
its  conduction  through  the  film  and  convection  along  the  rolling 
direction.  For  the  thermal  case,  the  program  computes  film  mean 
plane  and  surface  temperature  distributions  and  outputs  their 
average  and  maximum  values  over  the  contact. 

In  addition  to  traction  related  output,  the  program  computes 
the  contact  dimensions,  pressure,  deflection,  film  thickness, 
including  the  effects  of  starvation  and  inlet  heating,  Deborah 
number  and  a  surface  fatigue  index. 
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8 . 2  Program  Organization 

Computer  program  McFRIC  is  organized  into  9  modules  and 
comprises,  in  aggregate,  a  main  routine,  15  subroutines  and  17 
function  subprograms.  The  module  names  and  their  associated 
function  and  subroutine  names  are  listed  in  Table  8.1.  The  input 
data  comprises  34  elements  of  variable  data  and  17  run  control 
variables.  These  input  data  elements  are  contained  in  an  input 
file  prepared  and  named  by  the  user  by  editing  a  sample  file 
supplied  with  the  program  disk.  The  input  data  is  discussed 
further  in  Section  8.4.3  below. 

8 . 3  Program  Logic  Flow 

Figure  8.1  shows  the  overall  logic  flow  for  computer  program 
McFRIC.  On  execution,  the  program  requests  that  the  name  of  an 
input  data  file  be  typed  from  the  console.  This  file  is  then 
read  and  echoed  in  the  printed  output.  The  program  calculates 
the  contact  parameters  using  subroutine  HERTZ,  the  film  thickness 
(starved,  unstarved,  constriction  and  plateau)  using  subroutine 
FILM  and  Deborah  number  (subroutine  DEBORA).  The  microcontact 
parameters  are  computed  using  subroutine  GW  and  a  surface  fatigue 
index  is  computed  within  the  main  program.  The  program  enters  a 
loop  in  which  the  sliding  velocity  is  varied  over  a  user  spe¬ 
cified  range. 

For  each  value  of  the  sliding  velocity  the  fluid  traction 
forces  and  torque  are  computed  through  subroutine  TRACT  and  its 
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Table  8.1  MODULES,  FUNCTIONS  AND  SUBROUTINE  NAMES  FOR  MCFRIC 


MODULE 

MCFRIC 

INPOUT 

EHD 

GW 

TRACT 

STRESS 

FT  I  NT 

MAXAV 


FUNCTION 


SUBROUTINE 


INTEG,  SIMPS,  FI 


MCFRIC  (MAIN) 
PLOT 


F2 ,  F3 ,  FI 


OUTPUT 


OUTPUT 

INPUT 

PRINT 


SOLN  HERTZ 

ELINT 

FILM 


FUNO,  ERFC,  FUN 11  GW 

FUN321,  FUN1 ,  FUN32  DEBORA 

ANORM,  SPECEX ,  FUNC 
SPEC4 


TRACT 


FUN 


STRESS 

XSTRESS 


FTINT 


MAXAV 


HEAT 


HEAT 


subsidiary  subroutines.  TRACT  is  discussed  in  more  detail  below. 
At  the  completion  of  the  traction  computational  loop,  the 
remaining  program  output  is  printed.  At  the  user's  option 
printer  plots  may  be  specified  in  which  up  to  12  traction  related 
variables  may  be  plotted  against  sliding  speed. 

Figure  8.2  is  a  detailed  flow  diagram  for  subroutine  TRACT. 
As  noted,  it  computes  fluid  traction  point-by-point  over  the  con¬ 
tact  area,  at  a  specified  value  of  the  sliding  velocity  Au.  The 
contact  ellipse  is  divided  into  NY  strips  of  thickness  Ay  = 

2  b/NY.  The  index  ”1"  references  the  individual  strips.  For 
each  strip  the  half  width  a(y)  and  maximum  pressure  p(y)  are  com¬ 
puted  along  with  the  step  length  in  the  x  direction  Ax  =  2a(y)/NX 
(NX  =  number  of  elements  on  each  strip). 

The  index  "J"  references  each  of  the  NX  individual  rec¬ 
tangular  elements  along  a  strip,  in  the  x  direction.  The 

pressure  at  each  element  determines  the  viscosity  n(p)  and  the 

•  • 

limiting  shear  xc(p).  The  shear  rates  Yx  and  Yy  are  functions  of 
the  location  of  each  element,  i.e.,  the  x,y  coordinates  or, 
equivalently  the  I  and  J  indices.  The  stress  components  tx  and 
Ty  and  the  torque  contribution  of  the  element  are  then  found  as  a 
function  of  the  lubricant  properties  and  shear  rates  using 
subroutine  STRESS  which  implements  the  solution  scheme  outlined 
in  Section  6.0.  At  this  point,  for  an  isothermal  solution 
ITHERM  *  0,  TRACT  increments  the  index  J,  i.e., 
next  adjacent  element  along  the  Ith  strip. 


it  moves  to  the 


INITIALIZE  ARRAYS 
COMPUTE  ay 


1  1  0  CONTINUE 


COMPUTE  Fx.  Fy.  T. 
OVERALL  A  VC.  St  MAX. 
SOUD  St  FILM  TEMP. 


Fig.  8.2  -  LOGIC  DIAGRAM  FOR  SUBROUTINE  TRACT 
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For  a  thermal  solution  ( ITHERM= 1 ) ,  TRACT  calls  subroutine 
HEAT  (cf.  Fig.  8.2.1)  to  calculate  the  generated  heat  at  the  ele¬ 
ment  location  (I,J)  due  to  fluid  shear  (Of)  and  asperity  friction 
Qc.  An  initial  guess  of  the  maximum  film  temperature  rise  A6C  is 
made  and  given  the  temporary  variable  name  DTETA1 .  A  second  tem¬ 
porary  variable  DTETASl  is  set  to  zero.  The  surface  temperature 
at  location  J  is  now  estimated  as  outlined  in  Section  7.0  and 
then  set  equal  to  the  variable  DTETAS2.  Iteration  counters  ITERl 
and  ITER2  are  maintained  for  convergence  of  the  surface  and  film 
temperatures,  respectively.  They  have  limiting  values  MAXIT1  and 
MAXIT2  presently  set  within  the  program  to  50.  The  fluid  proper¬ 
ties  averaged  over  the  film  temperature  are  next  computed  as 
outlined  in  Section  7.0.  The  shear  stresses,  torque  contribution 
and  heat  are  then  recomputed  using  the  modified  properties.  A 
new  value  of  A0c(j)  is  then  computed  and  set  equal  to  the  tem¬ 
porary  variable  DTETA2 .  DTETA1  and  DTETA2  are  compared.  If  they 
differ  by  more  than  the  value  EPS1  (EPS1  =  0.001  and  EPS 2  =0.1 
are  the  values  presently  hard  coded  in  the  program),  DTETAl  is  set 
to  DTETA2  and  another  iteration  is  performed  for  A0c(x).  When 
convergence  is  achieved,  A6C(J)  is  set  to  the  last  value  of 
DTETA2.  Convergence  of  the  surface  temperature  is  then  checked 
by  comparing  DTETAS2  and  DTETASl.  If  it  has  not  occurred, 

DTETASl  is  set  to  DTETAS2  and  another  iteration  is  performed. 
Convergence  failure  messages  are  printed  to  the  screen  if  the 
maximum  number  of  iterations  is  exceeded  in  either  loop,  but, 


FROM  MAIN 
SUBROUTINE 


ITER1  -  ITER1  -t-  1 


DTETAS1  -  DTETXS2 


COMPUTE 
Q  -  Qo  OF 
ITER1  -  0 


COMPUTE 

»®<=0)  -  (Q»h-2)/8*K 
*«c(J)  -  DTETA1 


COMPUTE  A©s(j) 


TETAS2  - 


Fig  8.2.1  THERMAL  ANALYSIS 
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computation  proceeds  as  indicated.  (In  such  cases,  the  user  must 
then  exercise  discretion  in  using  the  results.  Convergence 
failure  has  not,  however,  been  a  problem  in  the  test  evaluations 
of  McFRIC  discussed  in  Section  10.0). 

When  convergence  is  achieved  for  each  element  of  a  given 
slice,  the  shear  stresses  and  differential  torque  are  integrated 
to  give  the  contributions  DFX(I),  DFY(I)  and  DT(I)  of  the  slice 
to  the  traction  force  components  FX  and  FY  and  to  the  torque  T. 

A  subroutine  MAXAV  is  also  used  to  find  the  average  and  maximum 
solid  and  film  temperature  at  the  slice.  When  all  slices  are 
complete,  an  integration  is  performed  over  the  slices  to  give  the 
resultant  values  of  FX,  FY  and  T.  MAXAV  is  then  reinvoked  to 
give  the  overall  average  and  maximum  values  of  the  solid  and  film 
temperatures. 

8. 4  User's  Information 

Section  8.4  is  a  detailed  description  of  how  to  use  the 
McFRIC  program.  Topics  in  this  section  include  instructions  for 
loading  and  executing  the  program,  and  a  discussion  of  the  output 
data  for  a  sample  case.  The  information  is  intended  to  be  suf¬ 
ficient  to  permit  a  user  to  successfully  run  and  interpret  the 
results  of  the  McFRIC  program. 

8.4.1  Installing  the  Program 

McFRIC  is  written  in  the  FORTRAN  language  and  is  executable 
on  IBM-PC-XT  or  AT  computers  using  the  DOS  operating  system.  It 
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is  recommended  that  the  program  be  executed  on  a  computer  which 
has  the  appropriate  Math  Co-processor  chip  (8087  or  80287) 
installed . 

The  program  is  delivered  on  two  floppy  diskettes.  The  first 
diskette  (Diskette  No.  1)  contains  the  following  files: 

MCFRIC.EXE  -  Executable  code  for  McFRIC  program 

MCF.BAT  -  Batch  file  for  executing  McFRIC  program 

SAMPLE . DAT  -  Sample  problem  data  file 

The  second  diskette  (Diskette  No.  2)  contains  source  code 
files  (.FOR)  for  the  modules  making  up  the  McFRIC  program: 

MCFRIC,  INPOUT,  EHD,  GW,  TRACT,  STRESS,  FTINT,  MAXAV,  and  HEAT. 

It  is  recommended  for  ease  and  speed  of  execution  that  the 
program  be  installed  on  the  IBM-PC  hard  disk.  To  copy  the 
program  onto  the  hard  disk,  use  the  following  procedure: 

1.  Turn  on  the  computer  and  set  the  default  drive  to  C: 

2.  Insert  Diskette  No.  1  into  Drive  A 

3-  Issue  the  command  COPY  A:*.* 

When  the  copy  operation  is  completed,  remove  the  floppy  diskette 
and  store  it  in  a  safe  location.  If  desired,  Diskette  No.  2  can 
be  copied  onto  the  hard  disk  using  the  same  procedure.  However, 
the  files  on  Diskette  No.  2  will  only  be  needed  if  the  program  is 


8.4.2  Executing  the  Program 

Once  the  program  has  been  copied  onto  the  hard  disk,  it  may 
be  executed  as  follows: 

1.  Set  default  drive  to  C  (C>) 

2.  Issue  the  following  command:  MCF 

The  program  should  begin  execution  by  displaying  a  heading  on  the 
screen  and  prompting  the  user  for  the  name  of  an  input  data  file. 
By  typing  SAMPLE . DAT ,  the  sample  case  data  file  can  be  used  to 
verify  that  the  program  is  executing  properly.  A  listing  of  the 
sample  case  data  file  is  given  as  Fig.  8.3.  The  corresponding 
program  output  is  included  as  Fig.  8.4.  Section  8.4.3  below  con¬ 
tains  instructions  for  creating  an  input  data  file  in  the  format 
required  by  the  program. 

8.4.3  Creating  an  Input  Data  File 

The  McFRIC  program  has  been  set  up  to  read  the  input  data 
from  a  standard  ASCII  file  which  can  reside  on  the  hard  disk 

(C: _ )  or  on  a  floppy  disk  (A: _  or  B: _ ).  The  recommended 

technique  for  creating  a  new  data  file  is  as  follows: 

1.  Use  the  DOS  COPY  command  to  copy  an  old  data  file  into  a 
new  file.  For  example,  by  issuing  the  command  COPY 
C: SAMPLE . DAT  C:CASE1.DAT,  a  new  file  named  CASE1.DAT 
will  be  created  on  the  hard  disk  containing  the  sample 
case  data. 
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THREE  TITLE  LINES  DESCRIBING  THE  RUN: 


TEST  CASE  CREATED  FROfl  HPAFB  TEST  DATA  FILE  M18N9gO 
LOAD  =  140.1  LBS,  ROLLING  VELOCITY  =  74B  IN/SEC 
AMBIENT  TEMPERATURE  =  183  DEG.  F 

MATERIAL  DATA: 


ELASTIC  MODULUS  OF  BODY  I  <lb/in#«2) 

[EMODUn  ..  29.5E6 

ELASTIC  MODULUS  OF  BODY  II  I3b/in*«2> 

[EMODCII  .,  29.5E6 

POISSON'S  RATIO  FOR  BODY  I 

I  PRAT  U)  I  ..  0.30 

POISSON'S  RATIO  FOR  BODY  II 

IPRAT !2! I  ..  0.30 

TENSILE  YIELD  STRESS  OF  SOFTER  MATERIAL  (lb/in#«2) 
I5I6MAEI  ...  3.62E5 

COULOMB  FRICTION  COEFFICIENT 


(EF)  . 

0.06 

DENSITY 

(lb/in»«3) 

'ROM!  . 

0.284 

SPECIFIC  HEAT 

(Btu/lba-deg.  F) 

ICM]  . 

0.11 

THERMAL  CONDUCTIVITY 

Ub/sec-deq.  F) 

t Akh]  .  7.34 

GEOMETRY  AND  ENVIRONMENTAL  DATA: 


RADIUS  OF  BODY  I  TRANSVERSE  TO  ROLLING  DIRECTION  (ini 
IRI11  .  0.71 

RADIUS  OF  BODY  I  IN  ROLLING  DIRECTION  (in) 

(RI2I  .  0.75 

RADIUS  OF  BODY  II  TRANSVERSE  TO  ROLLING  DIRECTION  (in) 
(Rill!  .  0.71 

RADIUS  OF  BODY  II  IN  ROLLING  DIRECTION  (in) 

[RII2T  .  0.75 

VELOCITY  OF  BODY  I  IN  ROLLING  DIRECTION  (ifl/secl 

tUIl  .  74B.18 

VELOCITY  OF  BODY  II  IN  R0LUN6  DIRECTION  (in/sec) 

tUill  .  748.18 

TRANSVERSE  DIRECTION  SLIDING  VELOCITY  (in/sec) 

tDELVl  .  0.0 

SPINNING  VELOCITY  (radians/sec) 

[ONEGA]  ....  0.0 

MENISCUS  DISTANCE  (in) 

tEJBJ  .  1.0 

CONTACT  NORMAL  LOAD  (lbs)  [PI 

.  140.1 

AMBIENT  TEMPERATURE  (dag.  F) 

ETAMB1  .  183.0 


Fig.  8.3  INPUT  FILE 


SAMPLE. DAI 
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LUBRICANT  DATA 


I 

i 

KINEMATIC  VISCOSITY  AT  0PERATIN6  TEMPERATURE  lest) 
r VISCSI  ....  1.31 
SPECIFIC  GRAVITY 
ISPGR]  .  1.7627 


TEMPERATURE  VISCOSITY  COEFFICIENT 
[BETA]  .  1.01E-2 

< 1/deg.  R) 

PRESSURE  VISCOSITY  INDEX 
[ ALPHA 1 3  ...  7. 14E-5 

(in**2/lb) 

THERMAL  CONDUCTIVITY 
[CONI  .  0.0086 

(lb/sec-deg.  F) 

SPECIFIC  HEAT 
[CFI  .  0.200 

(Btu/lba-deg.  F) 

AVERAGE  SHEAR  MODULUS 
[GEAR]  .  2443.69 

(lb/in««2) 

LIMITING  SHEAR  STRESS 
[TAUCI  .  12603.41 

(lb/in**21 

SURFACE  DATA: 


RMS  SURFACE  HEIGHT  FOR  BODY  1 

IRBII  .  2.0 

PUS  SURFACE  HEIGHT  FOR  BODY  II 

[R92I  .  2.0 

RMS  PROFILE  SLOPE  FOR  BODY  I 
CDELB1I  ....  0.03S 

RMS  PROFILE  SLOPE  FOR  BODY  II 
CDELB2I  ....  0.035 

LONER  FREQUENCY  LIMIT 

CFR11I  .  50.0 

UPPER  FREQUENCY  LIMIT 
CFR21I  . 10000.0 


(aicroinch) 

(aicroinch) 


ll/in) 

ll/in) 


RUN  CONTROL  DATA: 


INCREMENT  FOR  DELTAU  =  UI  -  UII  (in/sec) 

CUINCI  .  1.0 

MAI  I  MUM  VALUE  FOR  DELTAU  lilt /set) 

I  UMAX  1  .  10.0 

TYPE  OF  ANALYSIS  <0  *  ISOTHERMAL,  1  =  THERMAL) 

II THERM]  ...0 

PLOT  SELECTION  ARRAY  <0  =  NO  PLOT,  1  =  PLOT)  : 

PLOT  FLUID  SHEAR  FORCE  IN  R0LLIN6  DIRECTION  IFII? 
UVEC(l)]  ..1 

PLOT  FLUID  SHEAR  FORCE  IN  TRANSVERSE  DIRECTION  (FY)? 
(IVECC2)]  ..0 

PLOT  TORQUE  DUE  TO  FLUID  SHEAR  IT)? 

UVEC(I)]  ..0 

PLOT  EXPECTED  ASPERITY  TRACTION  FORCE 
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-  IN  ROLLINS  DIRECTION  <E(FX>)? 

[IVEC'4S 1  ..0 

-  IN  TRANSVERSE  DIRECTION  IF.(FY) )? 

CIVEC(5)J  ..0 

PLOT  EXPECTED  TORQUE  DUE  TO  ASPERITY  FRICTION  (E (T) ) 7 
t!VEC(6) I  ..0 
PLOT  TOTAL  TRACTION  FORCE 

-  IN  ROLLINS  DIRECTION  (FX+E(FX)>? 

r 1VEC (7)1  ..0 

-  IN  TRANSVERSE  DIRECTION  (FY+EIFY))? 

[IVEC(B)]  ..0 

PLOT  TOTAL  TRACTION  TORQUE  (T+E(T)P 
EIVECX9) I  ..0 

PLOT  C0I1B I  NED  FRICTIONAL  POKER  LOSS  (EL)? 

IIVECUO) I  .0 

PLOT  POKER  TRANSMITTED  THROUGH  CONTACT  (ET) ? 
tlVEC(tt)]  .0 
PLOT  LOSS  FACTOR  (LFI7 
EIVEC* 12) 1  .0 

NUHBER  OF  DIVISIONS  OF  CONTACT  ELLIPSE 

-  IN  ROLLING  DIRECTION 
CNXDIV I  ....11 

-  IN  TRANSVERSE  DIRECTION 
INYD1VI  ..,.41 


86 


R  ACT  ION  PREDICTION  **  AIR  FORCE  MATERIALS  LABORATORY  ♦»  TRACTION  PREDICTION 


KttltmtXmiltltltMllltUOXMKIIIH 

RCFRIC  COMPUTER  PROGRAM 
■  VERSION  :  1.0  RELEASE  DATE  :  B-15-Bfc 
-  SKF  TRJB0NET1CS  KING  OF  PRUSSIA,  PA 


TEST  CASE  CREATED  FROM  NPAfB  TEST  DATA  FILE  H18N9gO 
LOAD  =  140,1  LBS,  ROLLING  VELOCITY  =  748  IN/SEC 
A-BIENT  "EPiPERATURE  =  183  DEG.  F 


INPiM  CA'A  FILE  :  SAMPLE . DAT 


I.  H  A  T  E  R  I  A  L  VARIABLE 

E  NOD  Ml  -  YOUNG'S  NODUIUS  FOP  BODY  I  - 

EN0DI21  •  YOUNG'S  MODULUS  FOR  BODY  II  * 

PRAT ( 1 )  -  POISSON'S  RATIO  FOR  BODY  I 
PRAT (2)  -  POISSON'S  RATIO  FOR  BODY  II  = 
SIGNAE  -  TENSILE  YIELD  STRENGTH 

EF  •  COULOMB  FRICTION  COEFFICIENT  = 

RON  -  DENSITY 
CN  -  SPECIFIC  HEAT 
AKN  -  THERMAL  CONDUCTIVITY 


Fig.  8.4  McFRIC  OUTPUT  FOR 


S 

2.950000E+07  L8/IN»#2 
2.950000E+07  LB/IN«»2 
3.000000E-01 
3. OOOOO OE -01 
3.620000E+05  IB/IN«»2 
6. 000000E-02 
.2B4000E+00  LB/IN»»3 
. I 10000E+00  BTU/LBN-DEB.  F 
. 734000E+01  LB/SEC-DEG.  F 


FILE  SAMPLE . DAT 


H 


';A:‘T;os  FhESiC’Ol  **  «  I  R  r  0  R  C  E  MATERIALS  LABORATORY 


'ECT  CASE  CREATED  FROM  NPAFB  TEST  DATA  FILE  MlBNBgO 
i  'JD  t  140. 1  LOR,  ROLLING  VELOCITY  =  748  IN/SEC 
5KPVNT  TENPFRAT’JRE  =  183  DEE.  F 


I  .  8  E  3  1  E  T  R  Y  AND  ENVIRONMENTAL 

RU  •  RADIOS  OF  BODY  !  (TRANSVERSE  TO  ROLLING  DIRECTION! 

V.7  -  RADIOS  OF  BODY  I  (IN  ROLLING  DIRECTION) 

PHI  -  RADIUS  OF  BODY  II  : TRANSVERSE  TO  POLLING  DIRECTION) 
RL 12  -  RADIUS  OF  BODY  II  (IN  ROLLING  DIRECTION) 
ill  -  VELOCITY  OF  BODY  I  UN  POLLING  DIRECTION! 

U1I  -  VELOCITY  OF  BODY  II  UN  ROLLING  DIRECTION) 

DElV  -  TRANSVERSE  DIRECTION  SLIDING  VELOCITY 
OMEGA  -  SPINNING  VELOCITY 
EXB  -  MENISCUS  DISTANCE 
P  -  LOAD 

TAME  -  AMBIENT  TEMPERATURE 


TRACTION  PREDICTION  *t 


A  R  I  A  S  L  E  S 

=  7. lOOOOE-Ol  INCHES 
=  7.50000E-01  INCHES 
=  7.  lOOOOE-Ol  INCHES 
=  7.50000E-01  INCHES 
=  7.48180E+02  INCHES/SEC 
=  7.4B180E+02  INCHES/SEC 
*  .OOOOOE+On  INCHES/SEC 
-  .OOOOOE+OO  RADIANS/SEC 
=  t.OOOOOE+OO  INCHES 
=  1 . 40 100E+02  LBS 
=  .  18300E+D3  DEG.  F 
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m  trap- ion  predicidn  **  a  :  p  force  materials  laboratory  **  traction  prediction  »» 

T[ET  CAPE  CRPATEC  FROM  UFAFB  TEST  DATA  PILE  MlBNRgO 
1  "Ap  =  ;4D.I  LBS,  ROLLING  VELOCITY  =  748  IN/SEC 
ArAIENT  TF!*rERATLwE  =  1ST  DEG.  F 


III.  IP  BR  I  CANT  VARIABLES 

Vises  -  OIL  VISCOSITY  AT  OPERATING  TEMPERATURE 
SPGR  -  SPECIFIC  GRAVITY  OF  OIL 
FT AO  -  ABSOLUTE  VISCOSITY 
SETA  -  TEMPERATURE  VISCOSITY  COEFFICIENT 
ALPHA 1  -  PRESSURE  VISCOSITY  INDE* 

CCN  -  THERMAL  CONDUCTIVITY  OF  LUBRICANT 
GBAR  -  AVERAGE  SHEAR  MODULUS 

tau:  -  limiting  shear  stress 
cf  -  specific  hfat 


-  1.310000E+00  CENTI5T0KES 
=  1 . 762700E+00 

-  3.348248E-07  LB-SEC/IN**? 

=  t.010000E-02  1/DEB.  R 

=  7. 140000E-05  IN«2/LB 
=  8.600000E-03  LB/SEC-F 
=  2.443690E+03  LBS/IN**2 
=  I.2&0341E+04  LBS/IN«2 
=  . 200000E+00  9TU/LBH-DEG.  F 


^^t^ie^S^TRfimOI^ODEL  DEVELOP  HEMT(U)  NRCBEArTnG^KF^^^^^2/^^^ 
AEROSPACE  KING  OF  PRUSSIA  PA  J  I  HCCOOL  SEP  87 
AFUAL-TR-87-4079  FJJ815-84-C-5828 


UNCLASSIFIED 


F/G  11/8 


NL 


MICROCOPY  RESOLUTION  tEST  CHART 


*«  TR4CT  T3I4  FREDICT’ON  *•  AIR  FORCE  MATERIALS  LABORATORY  t«  TRACTION  PREDICTION  »# 


TEST  CASE  CREATED  FROM  NPAF6  TEST  DATA  FILE  H18N9gO 
LOAD  =  I  AO. 1  LBS,  ROLLING  VELOCITY  =  7A8  IN/ SEC 
AMBIENT  TEMPERATURE  =  183  DEG.  F 


IV.  SURFACE  ROUGHNESS  VARIABLES 

RO!  -  ROOT  MEAN  SQUARE  SURFACE  HEIGHT  FOR  BODY  l  = 

RO?  •  ROOT  HE  AN  SQUARE  SURFACE  HEIGHT  FOR  BODY  II  = 

RO  -  COMPOSITE  ROOT  MEAN  SQUARE  SURFACE  HEIGHT  = 
rr.L01  -  ROOT  MEAN  SQUARE  PROFILE  SLOPE  FOR  BODY  I 
DELQ2  -  ROOT  MEAN  SQUARE  PROFILE  SLOPE  FOR  BODY  II 
0F.L3  -  COMPOSITE  ROOT  MEAN  SQUARE  PROFILE  SLOPE  = 
'RII  -  LONER  FREQUENCY  LIMIT  * 

FR21  -  UPPER  FREQUENCY  LIMIT 


2.00000  MICROINCH 
2.00000  MICROINCH 
2. 82843  HICROINCH 
.03500 
.03500 
.04950 

50.00000  1/IN 
10000.00000  1/IN 
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»*  FACTION  PREDICTION  «»  AIR  FORCE  BATERIAIS  LABORATORY  •«  TRACTION  PREDICTION  •* 


TEST  CASE  CREATED  FROM  NPAFB  TEST  DATA  FILE  HIBNYgO 
LOAD  =  140.1  LBS,  ROLLING  VELOCITY  =  748  IN/SEC 
ABB I ENT  TEMPERATURE  =  183  DEB.  F 


V.  RUN  CONTROL  VARIABLES 

l1  INC  -  INCREMENT  FOR  DEITAU  =  1.00000  IN/SEC 
JKAK  -  MAXIMUM  VALUE  FOR  DELTAU  =  10.00000  IN/SEC 
[TWERfl  -  TYPE  OF  ANALYSIS  (0* I SOTHERRAL , 1 =THERRAL)  =  0 


VARIABLE  NAME  PLOT  SELECTION  INDEX 

U=PLQT  VS  U;  0=NQ  PLOT) 


FX 

1 

FY 

0 

T 

0 

EIFXI 

0 

E(FY) 

0 

EIT) 

0 

F>  ♦  E(FX) 

0 

FY  ♦  E(FY) 

0 

T  ♦  EIT) 

0 

EL 

0 

ET 

0 

LF 

0 

i 
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*«  TRACTION  PREDICTION  *»  AIR  FORCE  MATERIALS  LABORATORY  ««  TRACTION  PREDICTION  »• 


TEST  CASE  CREATED  FROM  NPAFB  TEST  DATA  FILE  M18N9gO 
LOAD  =  140. 1  LBS,  ROLLINS  VELOCITY  =  748  IN/SEC 
AMBIENT  TEMPERATURE  *  183  DES.  F 


I 

HERTZ  CONTACT  PARAMETERS 

AREA  -  AREA  OF  CONTACT  ELLIPSE 
DELTA  -  TOTAL  ELASTIC  APPROACH 
SI EM AO  -  MAT I HUM  HERTZ  CONTACT  PRESSURE 
A  -  CONTACT  ELLIPSE  SEMI AH  IS  IN  ROLLINS  DIRECTION  = 

B  -  CONTACT  ELL’PSE  SEH1AJIS  IN  TRANSVERSE  DIRECTION  = 


5.62472E-04  INCHES**2 
4.904B0E-04  INCHES 
3.73S19E+05  LBS/IN«»2 
1.36275E-02  INCHES 
1.11381E-02  INCHES 


«t  TRACI  ION  PREDICTION  »»  A  !  R  FORCE  H  A  I  F  R  1  A  1  S  l  A  B  0  R  A  T  0  R  Y  •*  TRACTION  PREDICTION  *» 


TEST  CASE  CREATED  FRQH  HP AFT  TEST  DATA  FILE  N18N9gO 
IPAD  =  140.1  LBS,  ROLLING  VFL0C1TY  =  74B  IN/SEC 
AMBIENT  'EHPERA7URE  1  183  DEG.  f 


FILM  THICKNESSES  AND  REDUCTION  FACTORS 


HCF  -  FULLY  FLOODED  CENTRAL  FILM  THICKNESS 
HHF  -  FULLY  FLOODED  CONSTRICTION  FILM  THICKNESS  = 
PHISC  -  CENTRAL  FILM  STARVATION  FACTOR  - 

HCS  ••  STARVED  CENTRAL  FILM  THICKNESS 
PHISM  -  CONSTRICTION  FILM  STARVATION  FACTOR 
HNS  -  STARVFD  MINIMUM  FILM  THICKNESS 
PHIT  -  INLET  HEATING  REDUCTION  FACTOR 
H  -  NET  CENTRAL  FILM  THICKNESS 


MB435E-06  INCHES 
3. B0854E-06  INCHES 
l.OOOOOE+OO 
6.48435E-0G  INCHES 
1 . OOOOOE+OO 
3.80854E-06  INCHES 
9.09233E-01 
5.B9592E-06  INCHES 


HH  -  NET  MINIMUM  FILM  THICKNESS 


S  3.46293E-06  INCHES 


»♦  FACTION  PREDICTION  »»  AIR  FORCE  MATERIALS  LABORATORY  H  TRACTION  PREDICTION  •• 


TEST  CASE  CREATED  FROM  NPAFB  TEST  OATA  FILE  f!18N9gO 
LOAD  =  140.1  LBS,  ROLLINS  VELOCITY  =  743  1N/SEC 
AMBIENT  TEMPERATURE  *  183  DCS.  F 


TRACTION  VARIABLES 

D  -  DEBORAH  NUMBER  'BASED  ON  =  .39B07E+03 

AVERA6E  PRESSURE  ) 

DELTA)!  -  LOAD  CENTER  OISPIACEMENT  =  .5302AE-04  INCHES 


««  TRACTION  PREDICTION  *»  AIR  FORCE  MATERIALS  LABORATORY  ••  TRACTION  PREDICTION  «• 


TEST  CASE  CREATES  FROM  NPAFB  TEST  DATA  FILE  MIBNRgO 
LOAD  =  140.1  LBS,  ROLLINS  VELOCITY  =  748  IN/SEC 
‘WSIFNT  TEWERATUPF.  =  183  OES.  f 


DEFINITION  OF  PARAMETERS 


DELTAU  *  SLIDING  VELOCITY  (IN/SEC) 

FI  -  FLUID  SHEAR  FORCE  IN  ROLLINS  DIRECTION  (LBS) 

FY  -  FLUID  SHEAR  FORCE  IN  TRANSVERSE  DIRECTION  (LBS) 

T  TORQUE  DUE  TO  FLUID  SHEAR  (IN-LBS) 

£FX  -  EXPECTED  ASPERITY  FRICTIONAL  TRACTION  IN  ROLLING  DIRECTION  (LBS) 

EFY  -  EXPECTED  ASPERITY  FRICTIONAL  TRACTION  IN  TRANSVERSE  DIRECTION  ILBS) 
ET  -  EXPECTED  TORQUE  DUE  TO  ASPERITY  FRICTION  (IN-LBS) 

EFUFX  -  TOTAL  TRACTION  FORCE  IN  ROLLING  DIRECTION  (LBS) 

EFYiFY  TOTAL  TRACTION  FORCE  IN  TRANSVERSE  DIRECTION  (LBS) 

EUT  -  TOTAL  TRACTION  TORQUE  (IN-LBS) 

EL  -  COMBINED  FRICTIONAL  PONER  LOSS  (HP) 

E1R  -  PONER  TRANSMITTED  THROUGH  CONTACT  (HP) 

LF  -  LOSS  FACTOR 


95 


ii  trait  c*  fRf Dir ' : 


i»  A  !  R  fORCE  MATERIALS  LABORATORY  »*  TRACTION  PREDICTION  u 


TEST  CASE  CREATED  FROM  NPAFB  TEST  DATA  FILE  «18N9gO 
:CAC  =  UC.l  LBS,  ROLLINS  VELOCITY  =  748  IN/SEC 
ft-BI ?«T  TEMPERATURE  =  183  OES.  F 


delta j 

M 

rr  T 

EFT 

EFY 

ET 

FXiEFX 

FYiEFY  TiET 

EL 

ETR 

LF 

.00^ 

.  O'EiOO 

.O'TE+OO  .OOE+OO 

•OOEiOO 

.00E*00 

.OOEiOO 

.OOEiOO 

.OOEiOO  .OOEiOO 

.OOEiOO 

.OOE+OO 

.OOE+OO 

1.000 

3.0IE+00 

.  TOE+0O  -1.60E-19 

1.86E-01 

•POE 400 

-9.96E-04 

3.19EI00 

. OOE+OO  -9.94E-04 

4.83E-04 

3.62E-01 

1 . 34E-03 

2.000 

5. 02E+00 

.  ■">'>£♦  00  3.34E-29 

1.04E-04 

. OOE+OO 

.OOEiOO 

5.02E+00 

.OOEiOO  3.34E-20 

1.52E-03 

5.69E-01 

2.L7E-03 

',.000 

5.98E+Cv 

. COE +00  -4.48E-19 

5.87E-0B 

.00E*00 

. OOEiOO 

5.9BE iOO 

.OOE+OO  -4.48E-19 

2.72E03 

6.7BE-CI 

4.01E-03 

*.000 

o.CE+OO 

. O'.? mTO  -J.51E-18 

3.30E-11 

. OOE+OO 

•OOEiOO 

fc.42E*00 

.OOE+OO  -1.51E-18 

3.B9E-03 

7.2BE-01 

3.35E-03 

5.000 

6.63'>0" 

. ’ TE+OO  -I.53E- 18 

1.86E-14 

.OOEiOO 

•OOEiOO 

4,63E<00 

.OOE+OO  -1.53E-18 

5.03E-03 

7.52E-01 

6.68E-03 

6.000 

e,75"iOO 

.OTFitV:  2.29E-18 

1.05E-17 

.OO'iOO 

.OOEiOO 

S.75E+00 

■OOE+OO  2.29E1B 

L.1JE-03 

7.65E-01 

B.02E-03 

\  TOC 

t>.61E+C0 

. OOE+OO  1.71E-18 

5.B8F-21 

. OOE+OO 

.OOEiOO 

L.81E+00 

•OOE+OO  1.71E-18 

7.22E-03 

7.72E-01 

9.36E-03 

9.  TOO 

S.3SEOO 

.OOEiOiT  -1.13E-18 

3.31E-24 

.OOEiOO 

.OOEiOO 

6.B5E+00 

.OOE+OO  -1.13E-1B 

8.  JOE-03 

7.77E-01 

1.07E-02 

9.000 

4.89E+00 

. OOC+CC  2.C3E-3B 

1.84E-27 

.00E+00 

.OOEiOO 

6.88EI00 

.OOE+OO  2.03E-18 

9.38E-03 

7.B0E-01 

1 . 20E-02 

10.000 

6.90E+00 

. OOE+OC  l.lfcE-18 

1.05E-70 

.OOF.  ♦  00 

•OOEiOO 

6.90E+00 

.OOEiOO  1.16E-1B 

1.04E-O2 

7.B2E-01 

1.34E-02 
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»*  TRACTION  PREDICTION  <»  AIR  P  0  R  C  E  MATERIALS  LABORATORY  **  TRACTION  PREDICTION  •• 

TEST  CASE  CREATED  FROM  NPAFB  TEST  DATA  FILE  «18N9gO 
LOAD  -  140.1  LBS,  ROLLING  VELOCITY  =  748  IN/SEC 
A“3!ENT  TEMPERATURE  =  183  DEG.  F 


LOAD  SHARING  DATA  AT  H  /  S  I  6  =  2.08452 

ZNCON  -  EXPECTED  NUMBER  OF  CONTACTS  =  .1877BE-25 

AVF  -  AVERAGE  CONTACT  FORCE  =  .64t87E-03  LBS 

TOTF  -  TOTAL  CONTACT  FORCE  .21555E+01  LBS 

AVA  -  AVERAGE  CONTACT  AREA  =  . U518E-0B  1N««2 

TOTA  -  TOTAL  CONTACT  AREA  =  .21627E-34  1N»«2 

DSUN1  -  NO.  OF  SUMMITS  ON  EQUIVALENT  SURFACE  =  .33215E+05 

N4  -  MEAN  SQUARE  CURVATURE  =  .47234E+07  IN«»(-2> 

R  -  ASPERITY  Tip  RADIUS  =  .30575E-03  IN 
K  -  SPECTRAL  EXPONENT  =  .10879E+0) 
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•<  TRACTION  PREDICTION  **  AIR  FORCE  MATERIALS  LABORATORY  •«  TRACTION  PREDICTION 


TEST  CASE  CREATED  FROM  NPAFB  TEST  DATA  FILE  HiBNfyO 
IPAD  --  140.1  LBS,  ROLLINS  VELOCITY  :  74B  1N/SEC 
AfBlE'.T  TEMPERATURE  -  183  DEB.  F 


SURFACE  FATIGUE  PARAMETERS 

ZNPLAS  -  EXPECTED  NO.  DF  PLASTIC  CONTACTS  =  .16354E+04 
APLAS  -  EXPECTED  AREA  DF  PLASTIC  CONTACTS  =  .499UE-09  IN**2 
PHIA  -  SURFACE  FATIGUE  INDEX  =  .81624E-06 


98 


»*  TRACTION  PREDICTION  *t  AIR  FORCE  MATERIALS'  LABORATORY  «•  TRACTION  PREDICTION  ** 


FI  VERSUS  OELTAU 


.OOOE+OO 
1. OOOE+OO 
2. OOOE+OO 
3. OOOE+OO 
4. OOOE+OO 
5. OOOE+OO 
6. OOOE+OO 
7. OOOE+OO 
8.000E+00 
9. OOOE+OO 
l.OOOE+Ol 


.OOOE+OO  fc.394E-01  1.379E+00  2.069E+00  2.758E+00  3.44BE+00  4.137E+00  4.827E+00  5.516E+00  4.206E+00  4.896E+00 


FI 


1 

1 

1 

1 

I 

1 
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2*  Use  a  commercial  editor  or  word  processing  package  (such 
as  WORDSTAR  in  non-document  mode  or  WORD  in  non- 
formatted  mode)  to  edit  the  data  in  the  new  file. 

As  already  noted,  the  required  input  data  file  structure  is 
illustrated  in  Fig.  8.3.  In  setting  up  the  file  structure, 
effort  has  been  taken  to  make  editing  of  the  input  data  possible 
without  reference  to  a  manual.  To  this  end,  all  of  the  headings, 
variable  descriptions,  units  and  symbols,  and  blank  or  dashed 
lines  shown  in  Fig.  8.3  are  built  right  into  the  data  file. 
(Important  Note?  Since  the  program  uses  fixed  format  FORTRAN 
read  statements,  it  is  vital  that  the  user  edit  only  the  data 
fields.  None  of  the  headings,  variable  descriptions,  blank 
lines,  etc.  should  be  deleted  or  changed  in  length  because  the 
program  expects  to  read  the  data  in  a  specific  location). 

The  input  parameters  are  organized  into  6  general  categories: 
title  information,  material  data,  geometry  and  environmental  data, 
lubricant  data,  surface  data,  and  run  control  data.  The  input 
data  fields  are  identified  by  the  characters  A  (Alphanumeric) 

F  (Floating  Point)  and  I  (Integer)  indicating  the  type  of  data 
expected.  Alphanumeric  data  consists  of  any  combination  of  num¬ 
bers  letters  or  symbols.  As  noted,  the  program  output  produced 
by  this  input  file  is  included  as  Fig.  8.4. 


9.0  TRACTION  TEST  PROGRAM 


As  part  of  the  contract  effort,  a  set  of  four  lubricants  were 
selected,  fluid  samples  were  acquired  and  a  traction  test  program 
was  designed  and  performed  within  AFML,  using  an  Air  Force  owned 
traction  rig.  The  dual  purpose  of  the  traction  test  program  was 
to  provide  input  data  for  McFRIC  and  to  evaluate  its  predictive 
accuracy. 

9. 1  Traction  Rig  and  Test  Variables 

The  parallel  axis,  two  disc  traction  apparatus  used  in  this 
study  is  shown  schematically  in  Fig.  9.1.  A  full  description  is 
given  in  Smith  [25] .  A  characteristic  of  the  rig  is  that  it 
requires  only  a  small  amount  (400  ml)  of  the  test  fluid.  The 
crowned  discs  are  separately  driven  via  independent  drive 
transmissions  and  variable  speed  motors.  A  normal  load  is 
applied  by  means  of  a  pneumatic  load  system.  A  variable  speed 
pump  is  used  to  circulate  the  oil.  A  heat  exchanger  in  the  inlet 
line  maintains  the  lubricant  at  the  desired  temperature. 

The  discs  are  made  of  52100  vacuum  degassed  alloy  steel  with 
a  hardness  of  60-63  Rc.  They  are  1.5  in.  in  diameter.  In  the 
test  program  the  disks  had  a  lateral  or  crown  radius  of  either 
0.71  in.,  providing  a  nearly  circular  contact  zone,  or  6.22  in. 
providing  a  4:1  aspect  ratio,  i.e.,  ratio  of  the  major-to-minor 
axes  of  the  contact  ellipse. 
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Fig.  9.1  AFML  TRACTION  RIG 


The  test  variables  employed  were: 

1.  normal  load 

2.  surface  speed  (rolling  velocity) 

3.  temperature 

4.  aspect  ratio  (1:1  or  4:1) 

When  a  test  is  performed  a  data  acquisition  system  records  the 
transmitted  torque  and  speeds  of  the  two  spindles.  The  raw  data 
is  subsequently  read  on  a  mainframe  computer  and  torque  is  con¬ 
verted  to  traction  coefficient,  and  the  rotational  speeds  are 
converted  to  sliding  speed  scaled  by  rolling  speed.  Fig.  9.2  is 
a  plot  and  Table  9.1  the  adjusted  numerical  data  from  a  typical 
test . 

9. 2  Lubricant  Selection 

The  lubricants  selected  after  joint  discussion  with  AFML 
were : 

1.  Mobil  RL714  base  stock,  a  synthetic  hydrocarbon 

2.  Monsanto  Santotrac  50,  a  traction  fluid 

3-  Halocarbon  Products  SAFETOL-R-3. 1,  CTFE  hydraulic  fluid 

4.  Montedison  Fomblin  Z,  a  perf luoroalkylether 

These  oils  were  chosen  for  their  chemical  diversity  con¬ 
sistent  with  being  of  intrinsic  interest  to  AFML. 

During  the  course  of  the  test  program  it  was  discovered  that 
Fluid  No.  4  had  an  adverse  effect  on  the  test  specimens.  Testing 
was,  therefore,  suspended  on  this  fluid. 
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TRACTION  DATA  FOR  MOBIL  RL-714  LUBE 

549  RPM  s  i_oad  =  55.5  lbs  (151392  pel)  >  T««pa80  F 
D I  as*e  t  e  r  *  1  . 500  In  Crown  Rad!l*b*3l  L  6-14  In 
TC (u)  Material  52100  Vac  Degas  Rc  60~63  Steal 
Data  File  K0808P1 


SLIP/ROLL 


Fig.  9.2  EXPERIMENTAL  TRACTION  CURVE 
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TABLE  9.1  ADJUSTED  TRACTION  DATA  POINTS 


TRACI  ION  DATA  FOR  HOB  It.  Rl-714  LUBF 
545  RHM  :  Load=55 .5  lb*.  (11,139?  pst>  :  Te«p=80 
D iamet er=l .500  In  :  Crown  Radii=6.3l  &  6.14  In 
Material  :  52100  Vac  Oegas  Rc  60-63  Steel 
Oata  File  K0808P1 


t 


SUP/ROLL 
♦.223247936558 
♦.214825449325 
♦.204092902270 
♦.189495552626 
♦.178735834729 
♦.170185391809 
♦.160154092165 
♦.149267595311 
♦.134189667323 
♦ • 1 2539b658076 
♦.118138554873 
♦.108107322673 
♦.091946167213 
♦.081469149176 
♦.0  71783689490 
♦.062609428466 
♦.047911851497 
♦.035990053183 
♦.026866994573 
♦.018988376614 
♦.006398940062 
♦.000000000000 
-.007104808654 
-.016589580465 
-.025369241611 
-.034220311411 
-.050478183774 
-.061473280307 
-.069725784006 
-.080001 743749 
-.094994514238 
-.105427538995 
-.114827477354 
-.124677192638 
-.135011059941 
-.148688057235 
-.158470937793 
-.167305547178 
-.1  7t>J30l  78568 
-.193437713521 
—  .2  0281 4  322Q  7  2 
-.21  1  (.46783243 


TRACTION  COFFF. 
♦.027674469013 
♦.027223613875 
♦.026952478819 
♦•026637741805 
♦•026205821547 
♦.025993459509 
♦•025806440712 
♦.025203151580 
♦.024448551508 
♦.023883652569 
♦.023285083033 
♦•022682200608 
♦•021420221116 
♦•020380410023 
♦•019564074996 
♦.018179835505 
♦•016020874931 
♦•013874833549 
♦•012077675182 
♦•009338437182 
♦•003345575069 
♦.000000000000 
-.003715537039 
-.008085418931 
-.010587427440 
-.012687189245 
-.015424554585 
-.016647504835 
-.01767*981813 
-.018523279248 
-.019916209915 
-.020679601353 
-.021301196610 
-.021769360058 
-.022631391080 
-.02368H925032 
-.02 392 6780125 
-.024241451783 
-.024387033208 
-  ■  0250  O'.'O  5226  2 
-.025277940731 
-.025753598940 
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9.3  Lubricant  Properties 


Shear  Modulus  G  and  Average  Limiting  Shear  Stress  tc 

The  extreme  conditions  of  pressure,  transit  time  and  shear 
rate  met  in  a  typical  EHL  contact  preclude  the  direct  measurement 
by  ordinary  laboratory  techniques  of  the  properties  G,  tc  and  the 
average  viscosity  n  under  relevant  conditions.  They  must,  there¬ 
fore,  be  deduced  directly  from  traction  measurements. 

Following  Tevaarwerk  and  Johnson  [4J  it  is  assumed  that  trac¬ 
tion  curves  of  the  type  shown  in  Fig.  9.3  are  available  at  or 
close  to  the  loads,  speeds,  aspect  ratios  and  temperatures  of 
interest.  These  tests  are  termed  simple  sliding  tests  because 
there  is  no  imposed  spinning  velocity  or  transverse  sliding  velo¬ 
city.  Only  the  sliding  velocity  in  the  rolling  direction,  Au,  is 
varied.  The  maximum  traction  coefficient  is  reached  at  high 
sliding  velocities.  At  sufficiently  high  sliding  speeds  the 
fluid  will  be  plastic  so  that  the  traction  force  Fx  will  be  the 
product  of  the  average  limiting  shear  force  Tc  and  the  contact 
area  *ab,  i.e., 

F x  58  *abTc  (9.1) 


Dividing  Fx  by  the  applied  load  P  gives  the  maximum  traction 
coefficient  umav: 


Using  Eq.  (9.2)  tc  may  be  computed  as: 

Tc  =  Mmaxp/1,ab  (9.3) 

For  high  Deborah  numbers  the  oil  will  act  purely  elastically 
at  low  sliding  rates.  Tevaarwerk  and  Johnson  [4]  show  that  in 
this  regime  the  resultant  shear  stress  integrated  over  the  ellip¬ 
tical  contact  area  gives  a  force  of 

Fx  =  8/3  (a2 b/Ph)  G  •  Au/u  (9.4) 

wherein  G,  the  composite  shear  modulus  of  the  fluid  and  con¬ 
tacting  surfaces,  is  assumed  to  be  constant  through  the  contact, 
u  =  Fx/P  is,  therefore,  linear  with  Au/u  and  has  slope: 

m  =  8/3  a2 b/Ph  •  G  (9.5) 

By  measuring  the  slope  ■  of  the  traction  curve  in  the  high 
Deborah  number  regime  one  may  thus  compute  G  as: 

G  =  3/8  Phm/a2 b  (9.6) 

Inasmuch  as  the  solid  surfaces  as  well  as  the  fluid  deform 
elastically,  this  elastic  modulus  value  represents  the  combined 
deformation  of  the  composite  fluid/solid  system. 

For  low  Deborah  numbers  the  oil  in  the  interface  behaves  as  a 
Newtonian  fluid,  i.e.,  the  shear  stress  is  related  to  sliding 
velocity  as: 
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t  *  n  Au/h 


(9.7) 


Considering  isothermal  conditions  but  using  Barus '  law  for 
pressure  dependence  of  viscosity,  Fx  is  obtained  as: 

F x  =  jj^dxdy  =  JJn(pjdxdy  •  Au/h  (9.8) 

=  *  a  b  n  Au/h 

=  u  u  a  b  TT/h  •  Au/u 

where  is  the  average  viscosity  across  the  contact  ellipse.  In 
this  case  u  =  Fx/P  also  be  linear  with  Au/u,  but  with  slope 

m  given  by 

m  =  u  *  a  b  n/hP  (9.10) 

For  intermediate  Deborah  numbers  the  traction  curve  at  low 
shear  rates  will  depend  jointly  on  the  viscosity  and  shear  modu¬ 
lus. 


At  intermediate  shear  rates  the  shape  of  the  traction  curve 
will  depend  jointly  on  the  values  of  n,  G  and  t~c,  with  their 
relative  influence  a  function  of  Deborah  number.  At  high  Deborah 
numbers  the  traction  curve  shape  at  intermediate  sliding  veloci¬ 
ties  depends  on  G  and  Tc.  At  low  Deborah  numbers,  intermediate 
sliding  velocity  behavior  depends  on  "n  and  t"c. 

One  conceivable  approach  to  estimating  the  lubricant  parame¬ 
ters  would  be  to  use  a  least  squares  fit  to  an  experimental  trac- 
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tion  curve.  In  this  approach/  n,  G  and  tc  are  estimated  as  the 
values  which  minimize  the  sum  of  the  squared  deviations  of  the 
predicted  and  actual  curves.  This  approach  could  be  prohibiti¬ 
vely  time  consuming  inasmuch  as  the  program  would  have  to  be  run 
repeatedly  for  various  trial  values  of  the  lubricant  parameters. 

As  a  practical  alternative  the  following  approach  was  adop¬ 
ted:  1)  t c  is  computed  from  Eq.  (9.3),  2)  G  is  computed  from 

Eq.  (9.5)  on  the  assumption  that  the  traction  curve  at  low 
sliding  speeds  is  elastic,  3)  n  is  taken  at  the  value  based  on 
its  ambient  viscosity  value  and  Barus'  equation  integrated  over 
the  contact  pressure. 

These  values  are  used  in  the  computation  of  the  Deborah 
number  as  discussed  in  Section  2.0.  If  the  Deborah  number  is 
large  enough,  n  is  irrelevant  since  the  initial  slope  depends  on 
G  only.  If  Deborah  number  is  small,  n  can  be  adjusted  as  needed 
to  agree  with  the  measured  slope  ■  using  Eq.  (9.10). 


9.4  Further  Lubricant  Properties 


In  addition  to  G  and  tc,  it  is  necessary,  both  for  the  design 
of  the  traction  tests  and  for  input  to  McFRIC,  to  know  six  pro¬ 
perties  of  each  lubricant  as  a  function  of  temperature.  These 
properties  are: 


viscosity  n 

temperature  viscosity  coefficient  tf 
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pressure  viscosity  coefficient  a 
specific  gravity  SPGR 
specific  heat  c 
thermal  conductivity  Kf 

A  computer  program  entitled  MATPROP  and  written  in  Fortran  77 
for  an  IBM  PC  and  compatibles,  was  developed  to  compute  the 
required  values  at  whatever  temperatures  T  (°F)  are  of  interest, 
for  the  three  oils  used  in  the  test  program.  The  relationships 
employed  within  MATPROP  are  given  below. 

9.4.1  Viscosity 

The  Walther-ASTM  equation  was  used  to  relate  viscosity  n  (cs) 
to  temperature  (°F)  as  follows: 


log(n+0.6)  =  exp[K-cinT] 


(9.11) 


where  K  and  c  are  fluid  dependent  constants.  By  substituting  n 
at  two  temperatures,  it  was  possible  to  solve  for  the  constants  K 
and  c  for  each  oil.  The  values  thus  determined  are  listed  below: 


»v 

r. 

B 

Lubricant 

K 

c 

1 

RL-714 

23.71805 

3.55330 

It 

Santotrac  50 

24.689883 

3.7017861 

CTFE 

43.66862 

6.82620 

| 

The  viscosity  values  used 

in  determining 

these  constants  were 

±  . 

taken  from  (26J,  1 3 7 J  and  128] 

respectively 

for  the  three  oils. 
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9.4.2  Temperature  Viscosity  Coefficient 

Given  the  viscosities  and  nQ,  corresponding  to  the  two 
temperatures  Tj  and  T0,  the  pressure  viscosity  coefficient  t*  is 
calculable  as 

«*  =  -^n(n1/n2)/(T1-T0)  (9.12) 

To  find  at  an  arbitrary  temperature  T,  Tj  was  taken  as 
T+20°  and  T0  as  T-20°.  and  n2  at  these  two  temperatures  were 

then  computed  from  Eq.  (9.11)  and  H  calculated  from  Eq.  (9.12)  as 

*  *  *n[n(T+20)/n(T-20) J/40  (9.13) 

9.4.3  Pressure  Viscosity  Coefficient 

RL714 

A  graph  of  u  vs.  T  supplied  in  [26),  was  fitted  by  the 
following  equation: 

a  =  exp  [6. 203294E-6  T2 -0. 00363T-8 . 89978 J  (9.14) 

SANTOTRAC  50 

Values  of  a  at  three  temperatures  were  given  in  [29]  for  a 
traction  fluid  known  to  be  Santotrac  50.  The  following  curve  fit 
was  developed: 


a  *  exp  [-0. 0057T  -  7.8858} 


(9.15) 


CTFE 


The  following  linear  approximation  was  developed  to  a  curve 
given  in  (31], 

a  =  1 . 82287E-4  -  6. 1598E-7*T  (9.16) 

This  equation  is  used  in  program  MATPROP.  It  does  not 

reflect  the  change  communicated  in  [30]  that  a  at  275°F  for  CTFE 

2 

is  4.276E-6  in  /lb.  The  program  is  believed  to  be  valid  for 
T<180°F. 


9.4.4  Specific  Gravity 
RL714 

The  equation  given  in  [26]  is: 

SPGR  =  8.177  -  1.13E-3  •  T  (9.17) 

The  following  equations  were  obtained  as  least  square  fits 
to  data  in  [27]  and  [28]. 

SANTOTRAC  50 

SPGR  =  0.9237  -  3. 45E-4T  (9.18) 

CTFE 


SPGR 


1.9303  -  9.3091E-4T 


(9.19) 


9.4.5  Conductivity 
RL714 


The  following  equation  was  given  in  [26]. 

K  =  ( 0. 813/ ( 12*SPGR) * ( 1 . 0-3E-4 (T-32)*0.2161  lb-sec/°F 
The  following  equations  were  obtained  as  least  squares 
SANTOTRAC  50 

K  =  0.061-1E-5T  ;  T  <  200°F 

=  0.065-3E-5T  ;  T  >  200 °F 


CTFE 

K  =  0.01016  -  81 . 647E-6T 

9.4.6  Specific  Heat 


RL714 

The  following  piecewise  linear  fit  was  obtained: 

c  =  1.175E-3T  +  0.405  T<200*F 

c  =  0. 210E-3T  +  0.598  T>200°F 

SANTOTRAC  50 

c  =  6.15E-4T  +  0.3840 


(9.20) 

fits: 


114 


CTFE 


The  value  c  =  0.20  was  used  independently  of  temperatures. 

9 . 5  Test  Design 

The  range  of  variables  used  in  the  tests  was  limited  by  the 
maximum  load  that  could  be  safely  applied  to  the  rig  (140  lbs), 
and  by  upper  and  lower  limits  on  the  rotational  speeds  (500  rpm 
and  10000  rpm).  Three  temperatures  were  chosen:  80°F,  180°F  and 
280°F.  Three  values  of  the  maximum  Hertz  stress  were  used:  50 
ksi,  206  and  366  ksi.  These  stresses  were  achieved  using  loads 
of  24.7  and  139  lbs  with  the  two  specimen  crown  radii.  Using  the 
lubricant  properties  relevant  to  the  temperature,  the  rolling 
speeds  were  calculated  using  computer  program  TRIBOS  [19]  so  as 
to  give  three  values  of  constant  film  thickness  for  each  tem¬ 
perature  regardless  of  load  or  aspect  ratio.  This  could  not 
always  be  achieved  because  of  the  physical  constraints  cited 
above . 

Tables  9.2  through  9.4  show  for  each  combination  of  load,  tem¬ 
perature  and  aspect  ratio  the  test  number  assigned  by  AFML,  the 
required  rolling  speed  (in/sec)  and  the  computed  film  thickness 
(win).  The  entire  set  of  tests  at  280°F  were  eliminated  for  the 
CTFE  fluid  because  of  concern  about  disk  damage  at  the  low  calcu¬ 
lated  film  thickness.  A  few  other  tests  had  to  be  discontinued  for 
various  reasons  accounting  for  some  missing  cells  in  the  test  matrices 
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10.0  TRACTION  DATA  REDUCTION.  ANALYSIS  AND  DISCUSSION 


10.1  Traction  Data  Base 

For  each  traction  test  performed  at  AFML,  the  value  of  Mmax 
was  read  from  the  printed  traction  data  output.  The  traction 
curve  slope  m  in  the  vicinity  of  zero  sliding  velocity  (Au-0) 
was  determined  using  the  method  of  least  squares  and  pairs  of 
values  of  Fx/P  and  Au/u  in  the  vicinity  of  Au=0. 

The  model  assumed  is: 

Yi  =  m*i+ei  (10.1) 

where  y^  is  the  traction  coefficient  measured  for  the  i-th  point, 
Xi  is  the  associated  slide-to-roll  ratio  (Au/u)i  and  is  a 
random  error  having  an  expected  value  of  zero. 

The  least  squares  estimate  of  m,  designated  m,  is  the  value 
which  minimizes  the  quantity: 

n  2 

s  =  l  (y i“mxi )  (10.2) 

i=l 

Setting  3s/3m  =  0  and  solving  gives, 

*  n  n  2 

m  «  l  xiY i/s  xi  (10.3) 

i=l  i =1 

A  program  called  TRSLOPE  was  written  in  BASIC  to  perform  this 
calculation  and  is  included  with  the  executable  Diskette  No.  1 
(cf.  Section  8.4.1). 


Table  10.1  is  a  data  form  used  to  compile  the  data  for  each 
traction  test  into  a  data  base.  With  the  data  base  software 
system  used  to  store  and  manipulate  these  data,  it  was  possible 
to  specify  certain  fields  by  a  formula  relating  the  value  of  the 
field  to  the  values  specified  in  one  or  more  of  the  other  fields. 
The  calculated  fields  are  denoted  by  an  asterisk  in  Table  10.1. 
The  formula  used  to  compute  these  fields  is  given  to  the  right  of 
the  colon  at  the  end  of  the  field  name.  Specifically,  the  shear 
strength  and  shear  modulus  are  computed  from  Eqs.  (9.3)  and  (9.6) 
under  the  supposition  that  the  fluid  behaves  elastically  in  the 
low  slip  region  and  plastically  at  high  slip.  It  is  further 
assumed  that  the  asperity  contribution  to  traction  is  negligible. 
Based  on  the  computed  film  thicknesses  relative  to  the  surface 
roughness  (roughly  2  uin)  this  was  almost  certainly  the  case  for 
all  but  a  few  of  the  CTFE  tests. 

Tables  10.2  through  10.4  give  the  complete  data  base  contents 
for  each  oil.  The  column  headings  correspond  to  the  field  names 
on  the  data  form,  but  have  been  truncated  in  some  cases. 

These  data  were  used  to  prepare  input  files  for  program 
McFRIC .  The  program  was  run  primarily  in  the  isothermal  mode 
unless  the  AFML  plots  showed  a  clear  indication  that  traction 
diminished  with  sliding  speed.  After  a  series  of  trial  runs 
showed  that  the  results  changed  negligibly  for  larger  values  of 
NX  and  NY,  the  number  of  divisions  of  the  contact  ellipse  were 
taken  as  NX  *  NY  *  41. 
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TABLE  10.1 


TRACTION  DATA  FORM 

LUBRICANT : 

FILE  NO. : 

ROTATIONAL  SPEED: 

*  ROLLING  VELOCITY.  0.0785  *  ROTATIONAL  SPEED 
LOAD: 

TEMPERATURE : 

*  ASPECT  RATIO:  SEMIAXIS  IN  ROLLING  DIR.  /  SEMIAXIS  TRANVERSE 
SEMIAXIS  IN  ROLLING  DIR. : 

SEMIAXIS  TRANVERSE: 

VISCOSITY: 

*  DYNAMIC  VISCOSITY:  VISCOSITY  *  SP  GR  *  0.0000001*5 
PRES  VISC  INDX: 

SP  GR: 

TEMP  VISC  INDX: 

CONDUCTIVITY : 

SP  HEAT: 

*  DENSITY:  SP  GR  *  62.5  /  1728 
FILM  THICKNESS: 

*  HERTZIAN  PRESSURE:  1.5  •  LOAD  /  (3.1*159  •  SEMIAXIS  IN  ROLLING  DIR.  •  SEMIAXIS  TRANVERSE) 

SLOPE  OF  TRACTION  CURVE: 

MAX  TRAC  COEF: 

*  SHEAR  MODULUS:  0.785*  *  FILM  THICKNESS  •  HERTZIAN  PRESSURE  •  SLOPE  OF  TRACTION  CURVE  /  SEMIAXIS  IN  ROLLING  DIR. 

*  SHEAR  STRENGTH:  HERTZIAN  PRESSURE  *  MAX  TRAC  COEF  /  1.5 
DEBORAH  NO: 

*  Computed  Automatically 
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TABLE  10.2  MOBIL  RL-714  DATABASE 


a 


s  i  aaaSSsHMHcnaaSSHHHHcc 

•  •  •••••••••••••••••••••••••«•• 

£  ■  ooooooooooooooooooocsooooooooo 

£  :  £3£33*SS8S8S8S32*223SSSSSSSS 

:! 331X332122128833138332222228 


0  0  «  a 


= : 3333383838338833S3833S8338838 
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TABLE  10.2  MOBIL  RL-714  DATABASE  (CONDT) 
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TABLE  10.3  SANTOTRAC -  50  DATABASE  (CONTD) 
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TABLE  10.4  HF-CTFE  DATABASE  (CONTD) 
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10.2  Predicted  Traction  Curves 

The  friction  force  Fx  predicted  by  MCFRIC  as  well  as  the 
measured  friction  force  were  plotted  against  the  sliding  speed 
Au,  for  all  of  the  experimental  traction  data.  The  complete  set 
of  these  curves  are  on  file  at  AFML. 

All  of  the  isothermal  fits  for  RL714  fluid  were  judged  to  be 
acceptable  even  though  some  of  the  tests  had  low  Deborah  numbers 
in  the  range  (9-22).  Figs.  10.1  and  10.2  are  examples  of  two 
distinctly  different  traction  curve  slopes  which  the  model 
nonetheless  appears  capable  of  accommodating.  Case  K18FOP1  had  a 
Deborah  number  of  12.14.  Case  K07BIV9  had  a  Deborah  number  of  466. 

Fig.  10.3  is  an  example  of  a  thermal  run  with  RL714.  It  is 
seen  that  the  program  grossly  overestimates  the  magnitude  of  the 
thermal  effect.  Fig.  10.4  shows  the  same  overestimate  of  the 
thermal  effect  for  a  CTFE  test.  The  program  was  accordingly 
modified  to  remove  the  effect  of  temperature  on  the  limiting 
shear  stress  while  leaving  the  effect  of  temperature  on  viscosity 
intact. 

The  revised  plots  are  shown  as  Figs.  10.5  and  10.6.  (The 
"bump"  in  the  predicted  curve  is  an  artifact  of  the  plotting 
program's  spline  fitting  routine).  The  fits  are  seen  to  be 
greatly  improved.  For  the  RL714  fluid,  there  remains  a  graphi¬ 
cally  discernible  thermal  effect  on  the  predicted  traction  curve. 
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The  experimental  traction  data  for  Santotrac  did  not  exhibit 
a  decrease  at  high  sliding  speeds  suggestive  of  a  thermal  effect 


Figs.  10.7  and  10.8  show,  respectively,  the  experimental  and 
isothermally  predicted  results  for  case  L28H7P1  which  had  the 
lowest  calculated  Deborah  number  among  the  Santotrac  tests 
(12.04),  and  for  case  L08B4W0  which  had  the  largest  Deborah 
number  value.  Both  fits  are  judged  to  be  acceptable. 

Fig.  10.9  shows  a  plot  of  the  results  of  a  thermal  run 
applied  to  Santotrac  case  L17H9P1  (Deborah  No.  =  179)  using  the 
original  version  of  the  program.  Fig.  10.10  shows  a  thermal  run 
with  the  modified  program.  The  program  output  does  show  a  tem¬ 
perature  decrease  at  high  sliding,  but  it  is  too  small  to  be 
graphically  perceptible. 

10.3  Viscoplastic  Regime 

The  most  serious  discrepancy  between  predicted  traction  and 
experimental  traction  force  values  occurred  in  the  9  CTFE  tests 
for  which  the  Deborah  number  was  less  than  unity. 

These  tests  were  all  performed  at  a  nominal  temperature  of 
180°F  at  which  the  pressure  viscosity  index  is  7.14E-5  in^/lb. 
Fig.  10.11  shows  the  isothermally  predicted  and  experimental  data 
for  File  No.  M18EOV7  for  which  the  Deborah  number  value  is  0.76. 
The  predicted  curve  appears  to  have  too  low  an  initial  slope  and 
fails  to  become  asymptotic  over  the  range  of  sliding  speeds 
represented  by  the  experimental  data. 
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As  noted  in  Section  9.4,  when  operating  in  the 


elastic/plastic  regime  (high  Deborah  Nos.)  the  model  essentially 
depends  on  the  two  values  (slope  and  maximum  coefficient)  that 
come  from  the  experimental  curve.  The  predicted  values,  there¬ 
fore,  cannot  help  but  reasonably  approximate  the  data.  For  low 
Deborah  numbers  the  predicted  curve  depends,  as  noted  in  Section 
9.4,  on  the  viscosity  n  averaged  over  the  contact  area  and  thus 
on  the  ambient  viscosity  and  more  problematically,  on  the  relation 
(Barus1  Law)  linking  pressure  and  viscosity.  An  approach  was 
accordingly  developed  to  correct  the  average  viscosity  to  yield 
the  empirically  observed  traction  behavior  in  the  low  Deborah 
number  regime  at  low  sliding  rates. 

Using  Eq.  (9.4),  one  may  calculate  n  giving  the  experimen¬ 
tally  determined  traction  curve  slope  m,  as: 

n  =  (h  P/u  n  a  b)  •  m  (10.4) 

Under  Barus'  Law, 


n 


where 


1/tt  ab 


a(y ) 

n(x,y)dxdy  = 


-b  -a (y ) 
p  r  a(y ) 

1/"  ab J  J  nQ  exp  [ap (x, y ) ] dxdy 
-b  -a (y ) 


P  (  x,  y  )  =  o0  [  1  -  ( x/  a ) 2  -  (y/b)Z]V2 


(10.5) 


(10.6) 


and 


a  (y )  =  a[l-(y/b)2]V2 


(10.7) 


Transforming  to 


Y  =  y/b 
and  X  =  x/a 

gives  n  =  nol}/*  (10.8) 


where 


1  +  [1-Y2]V2 


II  =  li(<*0o)  = 


// 


exp  [ao0[l-X2-Y2]l/2dxdy 


(10.9) 


-l  -[1-y2]1/2 


For  a  =  o,  Ii(o)  is  the  area  of  the  unit  circle,  i.e.  *,  so 
that  n  in  this  case  reduces  to 


II  was  evaluated  by  numerical  integration  over  the  range  of  aoQ 
from  1.0  to  100.  The  relationship  between  Ij  and  <*o0  was 
approximated  as  follows: 

Anlj  =  a  +  b  (aoQ)  +  c(aaQ)2  (10.10) 

The  values  of  a,  b  and  c  found  by  least  squares  over  four  subin- 
tervals  from  a aQ  *  1  to  50,  are  summarized  below: 
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RANGE  OF  ao0  £  b  c 


1-10  1.0681  0.7390  0.0099 

10-20  0.4215  0.8816  0.0018 

20-30  -0.0101  0.9257  7.150E-4 

40-50  -0.4848  0.9549  2.6146E-4 

Using  Eq.  (10.5)  in  (10.1),  the  value  of  nQ  required  to  make 

the  initial  traction  curve  slope  agree  with  the  measured  value  m 

is 


n0=  (hP/Ijuab)*m 


(10.11) 


Using  the  data  from  Table  10.3  for  File  M18EOV7,  the  product 
of  aoQ  is, 

aoQ  =  212,900  •  7.14E-5  =  15.2 

The  approximate  value  of  1^  may  be  computed  from  Eq.  (10.10): 
£nlx  =  0.4215  +  (.8816)  •  15.2  +  .0018(15. 2)2  =  14.24 
which  gives, 

I]L  =  el  4 . 24  =  1.525E6 

Using  this  value  of  1^  and  other  variable  values  from  Table 
10.3  in  Eq.  (10.11),  the  new  value  of  nQ  to  match  the  experimen¬ 
tal  traction  curve  slope  m  is 

„  _  3E-6  •  26  •  13.7 

no  “ 

1.52E6  •  252  •  7.78E-3  •  750E-3 


4.765E-8 


This  value  is  lower  than  the  value  3.35E-7  of  the  ambient 


viscosity  nQ  listed  in  Table  10.3  and  obtained  from  viscometer 
readings . 

The  average  viscosity  is  thus  not  the  source  of  the  low  ini¬ 
tial  slope  of  the  traction  curve  in  Fig.  10.11. 

10.4  The  Trachman-Cheng  Viscous  Model 

The  problem  appears  to  be  inherent  in  the  Trachman-Cheng 
nonlinear  viscous  law  which,  written  in  terms  of  the  average 
viscosity  and  yield  strength  is: 

t  =  n  /(  n/*c  +  i/7)  (10.12) 

—  • 

As  noted  t  approaches  tc  as  y  increases,  but  the  rate  of 
approach  appears  too  slow.  This  may  be  seen  in  Fig.  10.12  which 
is  a  replot  of  Fig.  10.11  with  a  wider  range  of  sliding  veloci- 
ties.  The  rate  of  change  of  t  with  respect  to  y  is,  in  fact, 

dx/dY  =  n  /  (n/*cY  +  l)2  (10.13) 

At  Y  =  0, 

dx/dY  =  n  (10.14) 

Even  if  the  slope  ‘n  is  correct  at  low  y  values,  it  decreases 
with  y  at  a  rate  which  can  be  quite  slow  if  n/xc  is  large. 

An  approach  that  might  be  usefully  pursued  is  to  develop  a 
lubricant  dependent  convergence  factor  k  such  that  the  above 
relationship  becomes 
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(10.15) 


X 


n/(n/tc 


+  1/kY) 


This  value  could  "tune”  the  approach  to  the  asymptotic 
limiting  stress  xc.  Such  a  factor  would  be  analogous  to  the 
parameter  x0  which  appears  in  the  Ree-Eyring  model  and  determines 
where  nonlinearity  sets  in.  It  is  recommended  that  this  or  com¬ 
parable  alternative  approaches  for  modulating  the  rate  of  con¬ 
vergence  to  xc  in  the  viscoplastic  regime  be  examined  in  future 
work . 

10.5  Alternative  Approach  to  Speeding  Convergence  Rate 

A  different  approach  that  does  not  involve  a  fundamental 
change  in  the  rheological  model  was  considered  for  restraining 
the  viscoplastically  predicted  traction  curve  from  becoming  asymp¬ 
totic  too  slowly  with  strain  rate. 

As  discussed,  the  empirically  measured  traction  curve  has  a 
slope  m  at  low  Au  and  approaches  an  asymptote  at  which  x  =  xc. 
Projecting  the  straight  line  portion  of  the  traction  curve  it  can 
be  deduced  that  the  traction  curve  becomes  horizontal  at  a  slide- 
to-roll  ratio  Au/u  in  the  vicinity  of: 

Au/u  =  umax/m  (10.16) 

The  shear  rate  at  this  point  is  thus 

Y  *  Au/h  *  u  wmax/mh  (10.17) 
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It  is  reasonable  to  require  that  t  given  by  the  T-C  relation 
has  risen  to  nearly  the  value  x  =  tc  at  this  shear  rate.  Taking 
|  x  =  0.99tc  and  using  Eq.  (10.17)  in  (10.12)  gives 

0.99  Yc  =  n/(n/Tc  +  rah/u  Mmax)  (10.18) 

Solving  for  n  gives, 

* 

n  =  99xc  m  h/u  ^maX  (10.19) 

Using  Mmax  =  Fx/P  where 

Fx  =  xc  •  rrab  (10.20) 

and  P  may  be  written 


P  =  2/3  oQ  •  itab  (10.21) 

so  that 

ymax  =  Tc/°o  (10.22) 

giving 

0  =  99  mhaQ  /(u  •  1.5)  (10.23) 

It  is  recalled  that  the  G  value  input  to  the  program  is 
deduced  from  the  empirically  measured  value  of  ■  under  the 
assumption  that  the  linear  portion  of  the  traction  curve  is  due 
»  to  elastic  deformation.  From  Gq.  (9.6),  m  may  be  re-expressed  as 

m  ■  8  a2  bG/3Ph  =  4  aG/*  h  oQ  (10.24) 


Using  (10.24)  in  (10.19)  results  in 

o  =  ( 8*99 )/3  •  (aG/»  u)  (10.25) 

Using  Eq.  10.5,  this  gives  the  modified  ambient  viscosity 
value  nG'of 

n0" =  264  aG/u  Ij 

Using  the  values  from  Table  10.3,  this  calculates  out  to, 

nQ'  =  264  —  hlMzUL- =  4.74E-6  (10.26) 

252  •  1.52E6 

Using  this  value  of  nQ  in  the  traction  calculation  (but  not 
in  the  film  thickness  calculation)  results  in  the  excellent  fit 
shown  in  Fig.  10.13.  A  recomputation  in  this  manner  has  been 
incorporated  into  McFRIC  and  is  used  wherever  the  Deborah  number 
falls  below  unity. 

The  nature  of  this  "fix"  is  revealed  by  Eq.  (10.25)  which  may 
be  written  as 

nu/aG  =  84 

This,  except  for  the  difference  in  computation  of  n,  is  the 
Deborah  number.  The  fix,  therefore,  amounts  to  changing  low 
Deborah  numbers  to  large  Deborah  numbers,  thereby  causing  the 
traction  curve  slope  at  low  slip  to  be  determined  by  the  elastic 
modulus. 
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Whether  this  is  correct  cannot  be  judged  with  the  present 
data.  It  is  necessary,  following  Johnson  and  Roberts  [14]  to 
measure  lateral  force  in  the  presence  of  spin  to  distinguish 
elastic  and  viscous  behavior.  As  it  stands  the  program  as 
modified  is  capable  of  predicting  the  results  of  simple  traction 
tests  for  all  three  test  oils. 

10.6  Estimating  G  and  tc 

McFRIC  requires  input  data  namely  G  and  tc,  that  must  be 
deduced  from  traction  tests  conducted  at  the  contact  conditions 
of  interest. 

To  overcome  this,  equations  were  fitted  to  the  traction  data 
using  stepwise  linear  regression  in  order  to  express  the  limiting 
shear  stress  and  shear  modulus  as  functions  of  temperature  T 
( °F ) ,  maximum  Hertz  pressure  p  (psi),  aspect  ratio  (AR),  and 
rolling  velocity  0  (in/sec). 

> 

\ 

j  The  coefficients  in  the  relationships  thus  determined  are 

i 

|  tabled  below  for  each  of  the  three  oils. 

I 

l 

i 

» 

» 

1 

i 

| 

! 
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LUBRICANT 

CONSTANT 

U  T 

AR 

E 

RL-714 

365.7 

1.05  -6.28 

379.  3 

0.0018 

Santotrac  50 

11576 

7.97  -53.4 

0 

0 

CTFE 

651.4 

-12.8 

1627 

0.0079 

Coefficients  for  G  Equation 

LUBRICANT 

CONSTANT 

U 

T 

AR 

E 

RL-714 

-1736 

-1.926 

-7.987 

0 

0.0334 

Santotrac  50 

3937 

-7.846 

0 

5538 

0.0322 

CTFE 

1599 

-7.017 

0 

6258 

0.0272 

Coefficients  for  tc  Equation 


Application  of  these  equations  to  the  data  used  in  their 
development  gave  reasonably  good  predictions.  If,  hovever,  G  or 
xc  was  small,  the  above  equations  occasionally  gave  negative 
values . 


The  equations  may  be  expected  to  be  most  reliable  for  con¬ 
ditions  which  fall  between  the  extremes  of  the  data  used  in  their 
development. 


« 
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